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PREFACE 


This  report  describes  a computer  program  designed  to  analyze  many 

•k 

aspects  of  the  geometric  performance  of  pseudoranging  navigation  satel- 
lite (navsat)  systems  for  users  either  on  earth  or  in  earth  orbit.  A 
navsat  system  includes  a fleet  of  satellites,  each  with  an  accurate 
clock,  which  transmits  emphemeris,  time,  and  other  signals.  These  sig- 
nals can  be  received  by  relatively  small,  inexpensive  equipment,  thus 
enabling  the  user  to  compute  his  position  and  time  accurately.  The 
NAVSTAR/Global  Positioning  System  (GPS)  currently  under  development  is 
such  a system.  The  overall  user  accuracy  of  such  a system  can  be  broken 
into  two  components  which,  when  multiplied  together,  yield  an  estimate 
of  the  user's  position  and/or  time  accuracy.  The  first  component, 
addressed  in  this  report,  depends  on  the  relative  geometry  among  the 
navsats  being  employed  and  the  user's  location.  The  second  Involves 
"system"  errors,  such  as  the  accuracy  of  the  ephemeris  data  of  the 
navsats,  propagation  effects,  clock  accuracies,  etc.  Convenient  com- 
puter analysis  of  the  geometric  performance  aspect  is  important  in  ad- 
dressing questions  of  alternative  orbital  configurations  for  the  navsats 
and  the  degradation  of  performance  due  to  failure  or  destruction  of  some 
or  many  of  them. 

The  impetus  for  this  research  came  from  Lieutenant  Colonel  Frank 
A.  Paparozzi,  Directorate  of  Space,  Hq  USAF  (AF/RDSA) , who  requested 
an  analysis  of  the  utility  and  feasibility  of  NAVSTAR/GPS  navigation 
support  to  high-altitude  satellites,  for  which  GPS  was  not  designed. 

The  research  • formed  as  a part  of  a Project  AIR  FORCE 

(formerly  Project  study  entitled  "Space  Warfare  Issues."  It 

should  be  of  use  to  military  and  civilian  defense  analysts  responsible 
for  the  design,  use,  and  survivability  of  GPS  and  other  U.S.  space- 
related  systems.  Additional  Project  AIR  FORCE  research  is  in  progress 
to  analyze  GPS  Phase  HI  design  features  needed  to  support  high-altitude 
space  navigation  and  user  equipment,  as  well  as  performance,  applica- 
tions, utility,  survivability,  and  alternatives. 

•k 

Ideally,  pseudoranging  is  a one-way  measurement  of  the  true  range 
plus  the  user's  unknown  time  oftset  with  respect  to  a master  reference. 
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SUMMARY 


This  report  presents  n discussion  (if  psciidoranging  navigation 
satellite  (navsat)  systems  (such  as  tlie  NAVSTAR/C I oba 1 Positioning 
System  (GPS)  satellite  system),  of  Geometric  Dilution  of  Precision 
(GDOP),  and  other  geometry-related  performance  parameters,  and  a com- 
puter progr.am  which  computes  them.  Included  are  satellite  selection 
algorithms  which  were  developed  to  minimize  the  computational  effort 
rec|uired  to  obtain  the  best  (or  nearly  best)  set  of  four  required 
satellites  for  computing  the  GDOP  for  either  a satellite-based  or  earth 
based  user. 

A computer  program  for  eartb-based  users,  developed  by  the  Aero- 
space Corporation,  was  acquired  and  extensively  modified.  Subject  to 
certain  constraints,  the  original  program  computed  the  number  of  navi- 
gation satellites  within  view  of  a user  at  any  location,  selected  the 
set  of  four  of  those  within  view  v/hich  would  minimize  navigation  errors 
and  computed  the  various  values  of  GDOP — namely,  the  three-dimensional 
position  error,  the  horizontal  position  i-rror,  the  altitude  error,  and 
the  time  error. 

At  Rand,  the  program  was  modified,  extended  to  accommodate  users 
in  any  earth  orbit,  and  optimized.  Further,  facilities  were  added  to 
give  the  user  convenient  and  powerful  input  and  output  control.  In 
short,  the  program  described  here  is  a flexible  "production"  program. 

in  addition  to  the  computation  of  navigation  satellite  coverage 
available  to  any  user  and  the  optimum  values  of  GDOP,  the  Rand  program 
also  includes  a feature  which  permits  the  variation  of  the  navsat 
antenna  beamwidth  and  determines  the  effect  of  this  variation  on  navi- 
gational accuracy  for  satellite  users.  There  are  no  restrictions  on 
the  shape  or  size  of  the  orbits  of  either  the  navigation  satellites  or 
the  user  satellite. 

The  computer  ()rogram  is  written  in  FORTRAN  IV  and  has  been  imple- 
mented on  an  IBM  170/1')8  computer  at  Rand.  Included  in  this  report 
are:  a program  listing,  an  ex|)ianatlon  of  the  variables,  a discussion 

of  the  operation  of  the  program,  and  sample  results. 
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' Pseudorangi ng  navigation  satellite  systems  transmit  one-way  signals 

I of  their  ephemeris,  time  (based  on  highly  accurate  and  stable  onboard 

I clocks  that  are  periodically  calibrated  by  ground  command),  and  other 

signals  that  allow  a user  with  fairly  simple,  lightweight  equipment  to 
I 

I accurately  determine  his  position,  velocity,  and  time.  Such  a system, 

1 called  the  NAVSTAR/Global  Positioning  System  (GPS),  is  being  developed 

by  the  military  services  at  the  Air  Force  Systems  Command  Space  and  Mis- 
sile Systems  Organization.  Tlie  GPS  will  eventually  consist  of  24  satel- 

* 

lltes  with  the  following  configuration:  three  orbit  planes  inclined 

63  deg  and  separated  by  120  deg  in  longitude;  eight  satellites  uni- 
formly distributed  in  each  plane  in  circular  orbits  at  about  10,900 
n ml  altitude  (12  hr  period).  (See  Fig.  1.)  The  system  is  being  de- 
signed to  provide  continuous  global  navigation  to  terrestrial  or  near- 

' earth  users  with  accuracies  on  the  order  of  tens  of  feet.  With  modlfl- 

I 

cations,  the  system  could  provide  high-altitude  satellites  (above  8000 
I n ml)  with  real  time  navigation  support.  This  navigation  support  could 

' effect  future  satellite  designs  and  operations,  tracking,  telemetry, 

j command  and  control,  and  many  other  space  applications. 

This  report  documents  a computer  program  which  simulates  the  orbital 
motion  of  a system  of  pseudoranging  navigation  satellites  (navsats)  and 

I ** 

the  motion  of  earth-based  users  or  a satellite-based  user.  The  program 
computes  the  number  of  navigational  satellites  within  view  of  the  user  and 
the  Geometric  Dilution  of  Precision  (GDOP)  values  which  are  dependent 
i primarily  on  the  user/satellite  geometry. 

^ The  position  accuracy  provided  by  such  systems  can  be  conveniently 

I divided  into  two  multlpllcaclve  factors — GDOP  and  other  "system"  errors. 

The  "system"  errors  depend  on  the  accuracy  of  the  ephemeris  data  and 
time  transmitted  by  the  navigation  satellites,  ionospheric  and  atmo- 
spheric effects,  and  various  mechanization,  electronic,  and  processing 

* 

The  GPS  Joint  Program  Office  is  considering  alternative  orbital 
conf igurat inns . 

** 

Ballistic  missiles  can  be  treated  as  satellites  with  orbital 
perigees  less  than  earth  Tadlus. 


i-rrors  in  tlif  navigation  satellite  ind  user  equipments.  Since  tiie  GDOP 


fac-Lors  depend  predominantly  on  the  user/navigation  satellite  geometries, 
they  can  he  analyzed  independently  of  system  errors,  which  depend  on 
cost,  technology,  and  effort.  This  allows  separate  analyses  of  alter- 
native orbital  configurations,  user  motion,  and  the  loss  of  some  of  the 
navigation  sateliites  by  interference  or  negation. 

The  original  i-ompnter  program  was  based  on  the  analyses  contained 
in  Ref.  1 and  was  designed  with  GPS  in  mind.  Both  the  referenced  report 
and  tile  original  computer  program  were  written  by  A.  Bogen  of  the  Aero- 
space Corporation.  The  original  program  was  designed  to  compute  the 
number  of  NAVS'FAR/GPS  satellites  available  to  an  earth-based  user  (or 
users)  and  the  various  GDOP  values  at  specified  times.  It  was  written 
so  that  the  set  of  navigational  satellites  can  be  changed  simply  by 
changing  the  initial  orbital  elements  of  each  satellite.  Satellite 
motion  is  assumed  to  be  Keplerian,  i.e.,  all  perturbing  forces  are  ne- 
glected. Both  satellite  coverage  and  GDOP  values  could  be  determined 
for  a single  earth-based  user  located  at  any  latitude  and  longitude  or 
on  a global  basis,  which  involved  an  assumed  uniform  distribution  of 
earth-based  users. 

At  Rand,  the  original  program  has  been  extensively  modified  and 
extended  to  permit  the  simulation  of  a satellite-based  user  and  any 
orbital  configuration  of  navsats.  There  are  no  restrictions  on  a satel- 


lite user's  orbit.  It  can  be  circular  or  elliptical  and  can  be  entirely 
above  or  below  the  altitude  of  the  navsats,  or  both  above  and  below 
(e.g.,  highly  elliptic  orbits). 

The  modifications  to  the  original  program  to  accommodate  a satel- 
lite user  include  algorithms  which  remove  unneeded  navigation  satellites 


that  are  in  view  of  the  satellite  user  in  order  to  decrease  the  amount 
* 

of  computation. 


A user  position  fix  requires  a determination  of  four  unknowns: 


three  components  of  position  plus  time,  thus  requiring  pseudoranging 

* A 

information  from  at  least  four  navigation  satellites.  Since  the  GDOP 


Some  combinations  of  satellite  user  orbit  and  navsat  orbit  configu- 
ration are  mori-  expensive  to  run  than  others;  this  will  be  explained  later. 


The  inclusion  of  an  explicit  time  calculation  reduces  the  .iccuracy 
and  stability  requirements  for  the  user's  clock,  thereby  making  t lu>  user 
equipment  simpler,  smaller,  lighter,  more  rugged,  and  cheaper. 
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I 


I 

I 


i 


I 


values  vary  [irlmarilv  with  the  relative  geometry  of  the  user  and  the 

navigation  satellites,  one  objective  is  to  select  that  combination  of 

four  navsats,  from  all  those  in  view  of  the  user,  which  will  yield  the 

minimum  (or  near  minimum)  values  of  GDOP. 

The  original  computer  program  cont. lined  an  algorithm  which  quickly 

selected  the  best  combination  of  four  (IPS  satellites  from  a maximum 

number  of  11  that  an  earth-based  user  could  see  with  a masking  angle 
* 

of  5 deg.  (See  Fig.  2.)  The  first  of  the  four  satellites  is  the  one 
nearest  the  user's  zenith  and  the  remaining  satellites  in  view  are  used 
three  at  a time  to  determine  the  best  combination  of  four.  This  is  an 
efficient  algorithm  that  results  in  the  selection  of  the  four  navsats 

'k'k 

which  will  yield  the  smallest  GDOP  almost  all  of  the  time.  For  11 
satellites  in  view,  the  maximum  that  GPS  would  provide  to  a terrestrial 
user,  120  combinations  of  three  are  possible  after  the  first  satellite 
has  been  selected,  and  there  are  330  different  combinations  of  four. 

For  a satellite-based  user,  assuming  sufficient  GPS  signal  strength 
in  all  directions,  the  masking  angle  depends  only  on  the  user's  altitude 
and  the  dimensions  assumed  for  the  earth  and  atmosphere.  Therefore,  as 
the  user's  altitude  increases,  so  will  the  number  of  satellites  in  view. 
For  a user  above  about  8000  n mi  altitude,  about  21  of  the  24  GPS  satel- 
lites would  be  in  view  if  the  satellites  were  uniformly  distributed. 

In  this  case,  after  selecting  the  first  satellite,  1330  different  combi- 
nations of  three  satellites  are  possible,  and  there  are  7315  different 
combinations  of  four. 

By  considering  the  geometry  of  the  satellite  user  relative  to  that 
of  the  satellites  in  view,  it  is  usually  possible  to  eliminate  about  half 


If  the  24  GPS  satellites  were  uniformly  distributed  on  the  spheri- 
cal surface  10,900  n mi  above  the  earth,  then  the  expected  number  in 
view  of  a user  with  a 5 deg  masking  angle  would  be  eight.  However, 
since  GPS  satellites  are  not  uniformly  distributed,  the  maximum  number 
in  view  can  be  11 . 

** 

Occasionally,  there  may  be  a navsat  with  a slightly  larger  zenith 
angle  that  would  result  in  a different  selection  of  the  other  three  nav- 
sats and  a smaller  GDOP.  A switch  is  provided  in  tlie  program  to  bypass 
the  algorithm  and  consider  all  combinations  of  four  navsats  in  view, 
but  is  far  more  expens iv<‘  to  use. 


navsat 


Local 

horizontal 


oi  tliem  from  further  consideration.  Thus,  tlie  original  program  was 
not  only  modified  t('  include  a satellite  user  hut  also  contains  an 
algorithm  for  such  a user  that  eliminates  navigation  satellites  in 
view  which  would  not  yield  small  (IDOP  values.  Versions  of  this  algo- 
rithm  could  be  used  in  satellite  user  equipment. 

In  general,  it  is  assumed  that  the  antenna  beamwidth  and  radiated 
power  of  the  navigation  satellite  are  adequate  to  assure  that  a satel- 
lite user  in  1 ine-of-sight  of  a navigation  satellite  can  receive  pseudo 
r.inging  information.  However,  tlu*  Rand  program  assumes  the  user  satel- 
lite has  an  omni-directional  antenna  but  permits  the  variation  of  the 
antenna  beamwidths  of  the  navsats  and  can  he  used  to  assess  the  result- 
ing effect  <in  satellite  user  coverage  and  CDOP  values. 

The  Rand  program  is  written  in  FORTRAN  IV,  and  is  implemented  on 
an  IBM  170/158  computer.  It  consists  of  a main  routine,  12  subroutines 
one  function  subprogram,  and  one  set  of  block  data.  An  external  sub- 
routine whicli  inverts  a symmetric  matrix  is  required.  There  are  three 
types  of  output:  Case  1 — user  on  a satellite;  Case  II — terrestrial 

user  .at  specified  Latitude  and  longitude;  and  Case  III  — terrestrial 
us(>rs  on  a net  of  Latitude  and  longitude  points  and  overall  global 
pi‘r  forma  nee. 

Section  II  discusses  the  CDOP  calculations  and  the  Rand  computer 
l>r(\’ram  in  gener.al  terms.  Section  III  describes  the  computer  program 
in  considerable  detail,  and  Section  IV  illustrates  the  more  convenient 
and  powerful  input  and  output  features  that  were  added  to  the  original 


progr.am. 


II.  DISCUSSION 


NAVICATIUN  SATELLITE  DEPLOYMENTS 

The  Rand  computer  program,  which  was  partly  written  with  GPS  Phase 
III  in  mind,  will  accommodate  up  to  a total  of  36  navsats  in  any  orbital 
arrangement.  One  additional  satellite  can  be  a user. 

Current  plans  for  the  fully  operational  GPS  Phase  III  call  for  24 
satellites  deployed  in  12-hr  (about  10,900  n mi  altitude)  circular  orbits 
inclined  63  deg  to  the  earth's  equatorial  plane.  The  24  satellites  will 
he  uniformly  distributed  in  three  orbit  planes  with  ascending  nodes 
separated  by  a longitude  of  120  deg.  (Alternative  orbital  configurations 
are  being  considered  by  the  GPS  Joint  Program  Office.)  The  eight  satel- 
lites in  each  plane  will  be  separated  by  45  deg,  and  the  phasing  of  the 
satellites  between  planes  will  be  chosen  to  provide  optimum  navigation 
geometry.  This  system  will  be  referred  to  frequently  in  the  following 
text  and  the  example  solutions  shown  later  are  based  on  this  system  of 
navigational  satellites. 

U S FJl  LOCATIONS 

Earth-based  users  can  be  located  at  any  earth  latitude  and  longi- 
tude. The  program  will  compute  the  GDOP  values  for  a single  user,  or 
for  multiple  users  on  a global  basis.  Currently,  the  program  will 
accommodate  only  one  satellite-borne  user.  There  aie  no  restrictions 
on  the  orbital  elements  of  any  of  the  navsats  or  the  satellite  user. 

NAVIGATION  SATELLITE  VISIBILITY 

In  Ref.  I,  it  is  stated  that  the  number  of  GPS  satellites  that 
would  be  visible  to  an  earth-based  user  would  never  be  less  than  6 
iior  greater  than  11  if  the  masking  angle  (elevation  angle  of  the  user's 
line  of  sight)  is  5 deg.  Larger  masking  angles  result  in  fewer  satel- 
lites being  visible  to  the  user  and  therefore  would  probably  decrease 
the  accuracy  in  determining  the  user's  position.  Although  smaller  mask- 
ing angles  would,  in  general,  increase  the  number  of  satellites  visible, 
atmospheric  effects  could  result  in  large  errors  in  the  pseudoranging 
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information  that  the  user  receives  from  GPS  satellites  near  the  user's 
horizon  due  to  propagation  errors. 

For  a satellite  user,  the  masking  angle  required  to  avoid  large 
atmospheric  propagation  effects  changes  rapidly  with  the  user's  alti- 
tude so  that  the  number  of  GPS  or  other  navigation  satellites  in  view 
increases.  In  fact,  a satellite  user  above  about  8000  n ml  altitude 
would,  on  the  average,  see  22  of  the  24  GPS  satellites.  At  this  alti- 
tude, only  two  satellites,  on  the  average,  would  be  obscured  from  the 
user's  view  by  the  earth  plus  200  n mi  of  atmosphere  above  the  earth. 
This  increase  in  the  number  of  GPS  satellites  visible  to  a satellite 
user  as  a function  of  the  user's  altitude  is  shown  in  Fig.  3. 


L 


USER  POSITION  ACCURACY* 


A user  of  a pseudoranging  navsat  system  needs  measurements  from 
four  satellites  in  order  to  determine  his  position  (three  components) 
and  time.  The  redundancy  afforded  by  having  more  than  four  satellites 
available  permits  a choice  of  the  one  set  of  four  which  will  yield  small 
values  of  GDOP.  These  values  depend  on  both  the  pseudorange  measurement 
errors  (system  errors)  and  the  relative  positions  of  the  four  GPS  satel- 
lites selected. 

The  geometric  relationships  between  the  user's  position,  the  posi- 
tions of  the  four  navigation  satellites,  and  the  four  pseudorange  mea- 
surements are  given  by  Ref.  1: 


where  1 

j 


b 


1,  2,  3,  4 and  is  the  index  of  each  of  the  four  equations, 
1,  2,  3 and  is  the  index  for  each  of  three  orthogonal  di- 
rections centered  at  the  user, 

the  j component  of  the  position  of  the  i*”^  navsat, 

the  component  of  the  position  of  the  user, 

the  pseudorange  measurement  from  the  user  to  the  l*"^  nav- 
sat, and 

the  user's  clock  bias  in  units  of  distance. 


These  formulations  are  partly  based  on  Information  from  Paul 
Jorgenson  of  the  Aerospace  Corporation. 
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Tlie  Information  noedod  to  solve  for  the  four  unknowns,  U.  and  b, 
are  the  pseudorange  measurements,  r.,  and  the  components  of  the  posi- 
tions of  the  four  navsats,  X...  If 

i.I 

= the  covariance  matrix  of  the  errors  in  the  pseudorange 
measurements,  and 

= the  covariance  matrix  of  the  resulting  errors  in  the 
navigation  parameters  (three  components  of  user  posi- 
tion and  user  clock  bias). 


then  it  can  be  shown  that 


C 

V 


C C C 


-T 


where  (I  is  a 4 x 4 matrix  of  tlie  coefficients  of  the  unknowns  in  the 

equations  obtained  by  linearizing  Eq . (1)  about  U and  b,  and  is  a 

function  of  only  the  direction  cosines  of  the  unit  vectors  from  the 

user  to  the  four  navsats  and  the  user's  c.lock  bias. 

There  are  some  small  geometric  effects  contained  in  C (for  ex- 

m 

ample-,  the  effect  of  elevation  angle  on  ionospheric  modeling  errors); 
these  effects  will  be  even  less  for  a satellite  user.  Thus,  a good 
apprt)x i ma t ion  of  the  geometric  performance  can  be  obtained  by  assuming 
that  C is  an  identity  matrix.  This  means  that  the  ranging  errors 

m on 

from  each  navsat  are  assumed  to  be  unity,  random,  and  uncorrelated. 
Hence,  to  a good  approximation, 

^ ■ (="'=)■'  • 

and  depends  only  on  the  direction  from  the  user  to  the  four  navsats 
and  is  in  no  way  dependent  on  the  distances  between  the  user  and  the 
satellites,  as  long  as  there  is  sufficient  signal  strength. 

The  various  GDOP  values  can  be  obtained  from  the  appropriate  ele- 
ments of  tlie  matrix  C^,  i.e.,  the  diagonal  elements.  If  Is  expressed 
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an  overall  measure  of  geometric  effect,  C1X)P,  is  obtained  from  the 
sc,  lare  root  of  t lu>  trace  of  the  matrix. 


CDOP  = 
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+ a 


22 
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44 
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where  j p ‘*22’  ‘^33  *■*’  errors  in  position  (x,  y,  z),  and  a^^  to 

the  user's  time  bias.  This  factor  includes  all  four  unknowns  (three 
dimensions  of  position  and  time)  and  is  the  conventional  measure  of 
overall  performance.  The  other  OOP  values  are: 


PDOP 

HIX)P 

V1X)P 

TDOP 

MDOP 


an  + a^2  ■•*33 


1/2 
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where  PIXIP,  HIXIP,  VDOP,  TDOP,  and  MDOP  are  the  multiplying  DOP  factors 
that  apply  for  the  three-dimensional  position  error,  the  horizontal 
position  error,  the  altitude  error,  the  time  error,  and  the  larger  com- 
ponent of  the  horizontal  position  error,  respectively. 

Because  the  effects  of  system  errors  (those  independent  of  geom- 
etry) are  multiplied  by  the  DOP  factors,  it  is  desirable  to  select  the 
four  navsats  of  those  available  that  will  yield  minimum  EX)P  values 
(generally,  PDOP). 

The  results  of  many  computer  runs  and  analytic  approximations  show 
that  there  is  almost  total  correlation  between  PDOP  and  the  volume  of 
the  tetrahedron  formed  by  lines  connecting  the  tips  of  the  four  unit 
vectors  from  the  user  toward  the  four  navsats.  Usually  (but  rarely 


The  origin  of  the  x,  y,  z coordinate  system  (a  right-handed  system) 
is  at  the  user's  position.  The  z axis  is  in  the  direction  of  the  user's 
vertical;  the  x axis  points  north  and  the  y axis  points  east  in  a plane 
normal  to  the  z direction. 

•irk 

For  example',  PDOP  can  be  shown  ti>  be  inversely  proportional  to 
six  times  the  volume  of  the  tetrahedron  and  directly  proportional  to 
the  trace'  of  a complicated  3x3  matrix. 


not),  the  larger  this  volume,  tlie  smaller  the  PDOP  values.  The  amount 
of  computational  time  required  to  compute  this  volume  is  much  less  than 
I the  computation  of  PDOP  itself,  which  involves  a matrix  inversion. 

' Thus,  the  program  was  designed  to  first  compute  the  volumes  of  the 

I tetrahedrons  associated  with  the  different  combinations  of  four  nav- 

] sats,  identify  the  combination  of  four  which  yields  the  largest  tetra- 

hedron volume,  and  then  use  that  combination  to  compute  the  DOP  values. 
Another  advantage  of  computing  the  tetrahedron  volume  as  a prelude  to 
computing  the  DOPs  is  that  the  time  rate  of  change  of  the  volume,  and 
consequently  an  estimate  of  the  time  rate  of  change  of  PDOP,  can  easily 
be  obtained; however,  this  is  not  implemented  in  the  program. 

SATELLITE  SELECTION 

An  earth-based  user  of  a 24-satellite  GPS  system  will  see  a maxi- 
mum of  11  of  the  24  satellites  if  his  masking  angle  is  as  low  as  3 deg. 

If  it  is  assumed  that  11  satellites  are  in  view,  and  the  one  closest  to 
tlie  user's  zenith  is  chosen  as  one  of  the  four  satellites  required,  there 
are  120  combinations  if  the  remaining  10  satellites  are  taken  three  at 
a time.  This  is  not  a large  number  of  combinations  to  be  investigated, 
so  the  original  program  did  not  include  an  algorithm  which  would  reduce 
the  number  of  combinations  to  be  examined.  There  is  an  option  in  the 
present  program  which  allows  the  satellites  to  be  combined  four  at  a 
time,  without  the  zenith  restriction,  but  still  using  the  tetrahedron 
vo 1 ume . 

In  the  case  of  a high-altitude  satellite  user  of  the  baseline  GPS 
configuration,  22  of  the  24  satellites  may  be  visible  at  any  time.  If 
I this  occurs,  there  are  1330  combinations  to  be  investigated,  and  there 

I is  an  option  in  the  program  which  allows  all  combinations  of  four  to 

I be  examined.  An  algorithm  was  develop'ed  to  eliminate  from  the  calcula- 

I tions  those  satellites  which  have  unfavorable  relative  locations.  This 

option  usually  cuts  the  execution  time  for  the  program  by  about  a factor 
of  10.  In  general,  about  half  of  the  satellites  in  view  can  be  elimi- 
nated. The  algorithm  is  derived  in  Appendix  A. 

Figure  4 is  a plot  of  PDOP  for  a satellite  user  in  a highly  elllp- 
I tic  orbit  versus  time  and  altitude,  and  the  variation  of  the  parameter 


Comparison  of  PDOPs  using  various  satellite  selection  criteria 


-14- 


valiie  for  C,  as  described  in  Appendix  A.  If  the  user  satellite  is  in  | 

a low  earth  orbit  (less  than  300  n mi)  or  above  the  CPS  altitude,  the  i 

1 

DOPs  will,  in  general,  be  the  same  for  all  of  the  satellite  selection  )i 

options.  However,  as  the  satellite  user's  altitude  approaches  the  al-  I 

titude  of  the  navsats  from  beneath,  the  zenith  selection  option  plus 

small  values  of  C will  not  yield  the  minimum  DOPs  because  the  algorithm  'i 

may  eliminate  some  navsats  which  should  be  retained.  This  is  caused  | 

by  tlie  fact  that  the  algorithm  retains  only  those  navsats  that  are  in  “ 

a band  determined  by  C which  is  normal  to  a line  from  the  earth's  center  | 

to  the  satellite  user,  rather  than  those  in  a band  normal  to  the  vector 

from  the  user  to  the  zenith  satellite.  However,  for  all  user  altitudes, 
the  use  of  all  navsats  in  view,  taken  four  at  a time,  will  yield  mini- 
mum DOPs  (based  on  using  the  tetrahedron  volume).  If  the  navsats'  or- 
bital configuration  involves  elliptic  orbits  or  more  than  one  circular 

orbit  "shell,"  there  are  cases  where  all  navsats  taken  four  at  a time  | 

must  be  used  or  extreme  care  exercised  in  selecting  a parameter  value  ) 

for  C.  See  Appendix  A for  details.  | 


HI.  COMPUTER  PROGRAM  OPERATION 


The  first  card  of  the  input  contains  the  switcli  for  ciuiosing  the 
case  desired  (I.OC:  l.OC  = 1,  satellite  user;  i.OC  = 2,  earth-based  user: 

LOG  = global  ca  I culittions)  ; the  number  of  navsats  (NJL)  ; and  the 
switch  for  choosing  the  navsat  selection  technique  (ISCMP:  ISCMP  = 

0,  zenith;  ISCMP  = I,  ail  satellites,  four  at  a time). 

There  are  three  types  of  calculations  performed  by  the  program: 
in  Case  1,  the  user  is  on  a satellite;  in  Case  II,  the  user  is  stationed 
at  a specific  latitude  and  longitude  on  the  earth's  surface;  and,  in 
Case  III,  a group  of  users  are  located  at  a set  of  latitudes  and  longi- 
tudes which  form  a net  over  the  whole  surface  or  a hemisphere  of  the 
earth . 

Along  with  each  of  the  three  cases,  there  is  a choice  of  the  nav- 
sat selection  technique  to  be  used.  In  Case  1,  when  the  user  is  on  a 
satellite,  one  choice  is  to  use  the  navsat  which  is  most  nearly  above  (>t‘ 
below  the  user  as  one  of  the  four  satellites  in  the  calculation  of  the 
volume  of  the  tetrahedron  (zenith  or  nadir  mode);  tlie  other  choice  is 
to  use  all  of  the  navsats  In  view,  taken  four  at  a time.  In  Cases  II 
and  III,  one  choice  Is  to  use  the  navsat  most  nearly  overhead  of  the 
earth-based  user  as  one  of  the  four  satellites  in  all  calculations  of 
the  volume,  and  the  other  choice  is  all  of  the  satellites  in  view  taken 
four  at  a time. 

In  all  cases,  the  basic  input  to  the  program  is  the  number  of  nav- 
sats and  their  orbital  elements:  eccentricity,  argument  of  perigee, 

right  ascension,  inclination,  initial  true  anomaly  (at  t = 0) , and 
period.  On  the  input  cards,  N.IE  is  the  number  of  satellites  in  the 
navsat  system.  P(N,K)  is  the  array  into  which  the  orbital  elements  are 
read,  where  N is  the  identification  number  of  the  navsat,  and  K = 1 to  S 
is  the  index  on  the  first  five  orbital  elements.  The  Index  K is  held 
c<instant  and  the  orbital  element  for  that  index  is  read  in  for  all  nav- 
sats (i.e.,  .ill  eci-en  t r i c i t ies  , then  all  arguments  of  perigee,  etc.). 

The  fin.il  inputs  for  the  n.ivsats  an-  their  periods,  which  are  read  into 
the  arr.iv  PER(N)  . After  tiu'se  inputs,  navs.it  user  related  data  is  read 
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For  Case  I,  NJL  + 1 is  the  number  of  the  user  satellite.  The 
first  satellite  user  data  card  of  a Case  1 input  contains  the  eccen- 
tricity, argument  of  perigee,  right  ascension,  inclination,  and  initial 
true  anomaly  of  the  user  satellite  at  t = 0.  The  second  card  contains 
the  period  of  the  user,  the  antenna  heamwidth  half-angle  of  the  navsats 
relative  to  a vector  from  the  navsat  to  the  center  of  the  earth  (AIN), 
and  a numerical  value  (C)  used  in  the  calculation  of  the  width  of  the 
band  on  the  spherical  surface  containing  the  navsats,  where  the  satel- 
lites for  consideration  in  the  "DOP"  calculations  will  be  sought  (see 
Appendix  A).  The  third  card  contains  the  time  increment,  in  minutes, 

at  which  the  calculations  are  desired  (INC),  and  the  total  number  of 

* 

these  time  increments,  plus  one  ( ITF) . For  each  run,  the  first  page 

of  output  contains  the  orbital  elements  and  all  other  input  values. 

The  subsequent  output  for  Case  I is  the  user  altitude,  VDOP,  hroP,  KDOP 

TDOP,  PDOP,  and  GDOP  at  the  chosen  time  increments,  plus  the  ilentifi- 

, cation  numbers  of  the  four  navsats  used  and  all  others  in  view.  If 

there  are  less  than  four  satellites  in  view,  there  will  be  no  print 

line  for  that  time  step — it  is  merely  skipped  and  the  program  continues 

The  first  navsat  user  data  card  for  a Case  II  run  contains  the 

latitude  (ATL) , longitude  (ONGL),  and  masking  angle  of  the  user  (ELEVAT 

see  Fig.  2).  The  second  card  contains  the  time  increment,  in  minutes, 

at  which  the  calculations  are  desired  (INC)  and  the  total  number  of 

* 

these  time  increments  required,  plus  one  (ITF).  The  output  for  Case 
II  is  identical  to  Case  1 — except  that  the  altitude  is  always  zero. 

In  Case  HI,  the  users  are  located  at  a set  of  latitudes  and 
longitudes  which  form  a net  over  the  surface  of  the  earth.  The  prime 
objective  of  this  Case  is  average  overall  system  performance  for  earth 
users  over  a period  of  time,  with  a secondary  objective  of  providing 
a "snapshot"  of  the  DOPs  at  specified  time  steps  at  each  longitude 
and  latitude  intersection  specified.  For  overall  performance,  global 
statistics  are  calculated  for  a uniform  distribution  of  users  over 
the  surface  of  the  earth  or  in  the  northern  hemisphere,  and  for  a full 


The  "plus  one"  accounts  for  the  "zero"  time  point. 

* ★ 

A uniform  distribution  of  users  is  approximated  by  the  DOPs  of 
users  at  a given  latitude  by  the  cosine  of  that  latitude. 
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orbital  period  or  a symmetrical  part  thereof.  The  total  time  input 
must  assure  this  condition  if  overall  performance  for  earth  users  is 
the  objective. 

The  net  of  users  can  be  chosen  in  various  combinations.  Basically, 
the  calculations  can  be  done  at  10  deg  steps  in  both  longitude  and  lati- 
tude, covering  the  whole  surface  or  the  earth;  or  at  10  deg  steps  of 
longitude  and  5 deg  steps  of  latitude  for  covering  the  northern  hemi- 
sphere. In  addition,  increments  of  these  basic  steps  can  be  chosen  so 
tliat  the  calculations  will  be  performed  at  every  20  or  every  30  deg, 
etc.,  of  longitude  and  at  every  20  or  every  30  deg.,  etc.,  or  at  every 
I or  every  15  deg,  etc.,  of  latitude.  Global  statistics  calculations 
averaged  over  time  are  output  at  the  completion  Cif  the  calculations,  but 
a "snapshot"  of  the  DOPs  at  each  longitude  and  latitude  selected  are 
available  at  any  time  point. 

The  first  navsat  user  card  for  a Case  III  input  contains  the  lati- 
tude step  size  (l.ATDEGl  and  the  masking  angle  of  the  users  (ELEVAT) . 

The  second  card  contains  the  latitude  increment  (LATIC)  and  the  longi- 
tude increment  (LONIC).  If  LATIC  = 3 and  LATDEG  = 5,  then  calculations 
will  be  made  at  0,  15,  30  deg,  etc.,  in  latitude.  If  LONIC  = 2,  then 
calculations  will  be  made  at  0,  20,  40,  deg.,  etc.,  in  longitude.  The 
second  card  also  contains  the  time  Increment,  in  minutes,  at  which  the 
calculations  are  desired  (INC)  and  the  total  number  of  these  time  incre- 
ments, plus  one,  required  to  complete  a full  orbital  period  or  to  reach 
a condition  of  symmetry  (ITF),  if  overall  performance  is  the  objective. 
The  final  two  variables  on  the  second  navsat  user  card  are  concerned 
with  the  "snapshot"  request.  If  IPFREQ  = 0,  no  "snapshot"  calculations 
are  made;  if  it  is  greater  than  0,  then  ITIME  is  set  to  an  Interger  and 
"snapshot"  calculations  will  be  made  and  printed.  For  example,  if  the 
user  has  cliosen  INC  = 10  and  ITIME  = 5,  then  "snapshot"  output  will 
occur  at  time  = 0,  50,  100,  etc.  If  the  "snapshot"  calculations  are 
requested,  the  user  will  always  get  output  at  time  = 0. 

Since  the  program  is  quite  costly  for  global  calculations,  the 
computation  timi'  step  should  be  small  enough  to  im'wnc  a representative 
global  distribution.  The  sizi>  of  the  time  step  will  depend  on  the 
navsat  orbital  configuration.  Output  for  the  av«'rage  overall  global 
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per f ormiuice  are  elev;ition  distribution,  latitude  elevation  distribution, 

I aecumulative  global  distribution  of  DOPs,  maximum  and  minimum  number 

of  satellites  seen  at  tlu'  latitudes  and  longitudes  of  the  net,  the  prob- 
j ability  of  seeing  exactly  'A  satellites,  and  the  probability  of  seeing 

■ N or  more  satellites  both  by  latitude  and  on  an  average  (weighted  by 

I the  cosine  of  latitude)  global  basis. 

Table  I (pages  19-22)  contains  an  explanation  of  the  variables  in 

■ the  MAIN  program,  and  Table  2 (pages  21-27)  contains  an  explanation  of 
the  variables  in  the  subroutines  and  in  their  calling  sequences.  Table 
'3  (pages  28-29)  shows  the  manner  in  which  the  inpu.  cards  are  .set  up. 

Appendix  B provides  a listing  of  the  entire  program. 

A few  timed  runs  were  made.  Table  4 (page  30)  shows  a comparison 
of  some  of  the  computing  parameters. 

The  program  has  been  dimensioned  to  accommodate  a total  of  36  nav- 
sats.  In  Case  I,  this  means  36  navsats  and  one  satellite  user.  Cases 
11  and  111  are  limited  to  a maximum  of  36  navsats.  All  cases  can  evoke 
the  option  of  using  all  navsats  in  view,  taken  four  at  a time,  but,  in 
i Case  in  and  in  Case  1 (for  high-altitude  satellite  users),  the  run  | 

j time  would  be  prohibitive.  If  less  than  36  satellites  are  to  be  used, 

the  dimensions  of  matrix  KXX(1,4)  (in  MAIN,  and  in  subroutines  TAAT  ‘ ' 

and  TAl.l.)  should  be  reduced  to  cut  down  core  storage.  The  index  T is  i 

I 

calculated  from  the  formula  i 

r!  (n  - r)!  ’ 

where  n = number  of  navsats  and  r = 4. 

For  Case  I,  t lu'  calculations  become  less  accurate  as  the  altitude 
of  a user  above  the  navsats  increases.  The  subroutine  which  inverts 
the  OOP  matrix  is  in  double  precision  in  order  to  alleviate  some  of 
this  problem.  However,  above  a ratio  of  satellite  user  altitude  to 
navsat  altitude  of  about  20,  the  matrix  becomes  sufficiently  ill- 
conditioned  that  errors  haejin  to  occur  in  the  DoPs. 

I I 


I i 

1 


EXPLANATION  OF  VARIABLES  IN  MAIN  PROGRAM 


Table  1 (continued) 


Table  1 (continued) 


This  portion  of  the  program  is  explained  in  Ref. 


Table  1 (continued) 


I 
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Tabie  2 

EXPLANATION  OF  VARIABLES  AND  CALLING  SEQUENCES  OF  SUBROUTINES 


The  subroutine  call  with  the  variable  names  in  its  argument  list,  as 
it  appears  in  MAIN,  is  shown  first.  The  subroutine  name,  with  the 
variable  names  in  its  argument  list  and  an  accompanying  explanation, 
appears  second. 


ORBINI(N) 


ORBINI (1) 

P(I,2)=P(I,2)*C(1) 
P(I,3)=P(I,3)*C(1) 
P(I,4)=P(I,A)*C(1) 
P(I,5)=P(I,5)*C(1) 
P(I,9)=sln  i 
P(I,10)=cos  i 
P(I,ll)=sin  n 


I,  satellite  number 
Argument  of  periapsis  x 0.017A5 
Right  ascension  of  the  node  x 0.01745( 
Inclination  x 0.01745 
Initial  true  anomaly  x 0.01745 
1,  Inclination 

12,  right  ascension  of  node 


P(I,12)=cos  12 

P(I,21)=a,  orbit  semi-major  axis  in  feet 
SAA=earth  synchronous  satellite  radius  in  feet 
P(I,23)=ratio  of  synchronous  radius  to  orbit  semi-major  axis 
P( I , 6)=a( l-e2)  e,  eccentricity 


P(I,7)=  ^l,  gravitational  constant 

P(I,8)=  Vv/[a2(l-e2)] 


Calculated  only 
to  change  units 
from  degrees 
to  radians 


POINT(LONG,LAT,TID,UPV) 

POINT(ALO,ALA,TIM,VEC) 


ALO,  longitude 
ALA,  latitude 
TIM,  time  in  days 

VEC,  vector  from  the  center  of  che  earth  to 
a user  at  a specific  latitude  and 
longitude 


EW=long(rad)+27T(  t ime) 
SN=lat (rad) 


I 


Table  2 (continued) 


VEC(1)= (radius 

of 

the 

earth 

at 

the 

equator) (cos 

lat) (cos 

long) 

VEC( 2) = (radius 

of 

the 

earth 

at 

the 

equator) (cos 

lat ) (sin 

long) 

VEC(3)=(radius 

of 

the 

earth 

at 

the 

equator) (sin 

lat) 

COVNAV(G,NSGD,4  SIGT) 

COVNAV(G, ID,NAT,SIG)  G,  direction  cosine  matrix 

ID,  index  of  the  navsats  which  were  chosen 
as  the  "best"  four 

NAT,  tile  number  of  them  (always  four) 

STG,  dilution  of  precision  parameters  re- 
turned to  M\1N 

This  subroutine  is  the  only  one  wliicti  requires  a local  system 

T -1 

subroutine — a matrix  inversion.  This  subroutine  calculates  [G  G] 
and  returns  the  values  in  SlG(l-6),  whicli  are  the  DOPs . 


TRMATX(TRANS,I) 

TRMATX(TR,1) 

TR(  1 ,1  ) = cos  f!  cos 
TR(l,2)=-cos  n sin 
TR(1  , 3)=sin  SI  sin 
TR(2,l)=sln  -Q  cos 
TR(2,2)=-sin  sin 
TR(2,3)=-cos  sin 
TR( 3, l)=sin  1 sin 
TR(3,2)=sln  i cos 
TR(3,3)=cos  i 


TR,  3x3  coordinate 
1,  satellite  number 
u - sin  U cos  i sin  u 
u - sin  fi  cos  i cos  u 
i 

u + cos  SI  cos  i sin  u 
u + cos  SI  cos  1 cos  u 
1 
u 
u 


transformation  matrix 


u,  argument  of 
latitude 


< 

; 

i 


MATMUL(TRANS,Q,R) 


t 
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I'able  2 (continued) 


MATMUL( TRANS. QC, AC) 

MATMI)L(T,V,0)  T,  3 x 3 coordinate  transformation  matrix 

V,  vector  to  be  transformed 
0,  vector  returned 

VOLUME(USUV,ISIC,VOLl)M) 

VOLUMfUUVRC.IDSAT.VOL)  UVEC,  unit  vector  from  the  user  to  each 

navsat  In  view 
IDSAT,  Index  of  the  navsat 
VOL,  volume  of  the  tetrahedron  formed  by 
the  set  of  four  satellites  used 

0RB1T(N,TID.PER,RF,VE,AA) 

ORBIT(I,T,PER,R,VEL,AC)  N,  I,  satellite  number 

TID,  T,  time  in  days 
PER,  period  of  a satellite 

RE,  R,  vector  to  a satellite  from  the  center 
of  the  earth 
VE,  VEL,  velocity 

AA,  AC,  acceleration  due  to  gravity 

This  subroutine  iterates  for  the  eccentric  anomaly  and  computes 
the  true  anomaly.  If  time  = 0,  the  following  is  computed: 

V(I)=0P(I,5)=v,  true  anomaly 
SINE=  [ ^1-e^  sin  v]/(l  -f  e cos  v],  sin  E 
E(I)=E,  arcsin  E,  eccentric  anomaly 
BIGT(I)=time  from  perigee  for  satellite 

If  time  > 0,  E(I)  is  obtained  by  iterating,  then  the  following  is 
computed : 

SINV=  (yl-e^  sin  E]/[l  - e cos  E)  = sin  v 
V(I)=arcsin  v,  true  anomaly 

U( I )*V( I )-K)P( 1 , 2) , V + argument  of  per laps  Is  = v + w * u,  argu- 


!not  used 
at  present 


ment  of  latitude 


Table  2 (continued) 


0P(I,13) 
0P(I,14) 
0P(I,15) 
0P(1 , 18) 
Q(l)  = I 
Q(2)  = 0 
Q(3)  = 0 
QEL(l)  = 
gEL(2)  = 
QEI.O)  = 
QC(1)  = 
Qt;(2)  = 
QC(3)  = 
OP(l,24) 
OP(l , lb) 


= cos  u 
= sin  u 
= latitude 

= u,  argument  of  latitude 
a(l-e^)/[l  + e cos  v] 


sin  V j 

U /'JaXT^e^) 


0 


(1  + e cos  v)^  u/(a^(l-e^)‘ 


0 

0 


= ear  til  rotation 
= longitude  of  the  satellite 


for  calculating  R 


for  calculating  VEL 


for  calculating  AC 


TAAT(NSS,KCOM, KXX) 

TAAT (MAX, MXX , MATRIX)  MAX,  1 less  than  the  number  of  navsats  to 

be  examined  for  use 

MXX,  the  number  of  combinations  of  those 
navsats 

MATRIX,  contains  indices  of  the  navsats  to 
be  examined,  three  at  a time 

This  subroutine  sets  up  the  sequence  of  navsats  which  are  to  be 
examined,  using  either  the  one  above  or  below  the  user  as  one  of  the 
four  in  eacli  calculation  of  the  tetrahedron  volume. 


TALL (KOT,KCOM, KXX) 

TALL(MAX, MXX, MATRIX)  MAX,  the  number  of  navsats  to  be  examined 

for  use 

MXX,  the  number  of  combinations  of  those 


1 


navsats 
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Table  2 (continued) 

MATRIX,  contains  Indices  of  the  navsats 
to  be  examined,  four  at  a time 

This  subroutine  sets  up  the  sequence  of  navsats  which  are  to  be 
examined,  using  all  satellites  taken  four  at  a time. 

B IU)CK_  DATA 

Contains  various  parameter  values  used  in  the  program,  although 
not  all  of  them  are  used  In  the  current  version. 

VECT0R(V1,I,V2,V3) 

This  subroutine  performs  vector  additions,  subtractions,  and 
cross  products. 

UNIVEC(U,UV) 

This  subroutine  calculates  unit  vectors. 

D0T(V1,V2) 

This  function  subprogram  calculates  dot  products. 
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Table  3 (continued) 


The  basic  'Lcngitude  step  is  always  10  deg.  LATDEG  = 10  is  for  northern  and  southern  hemisphere 
coverage;  LATDEG  = 5 is  for  northern  hemisphere  coverage. 
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Table  4 

COMPUTING  COSTS  (IBM  370/158  COMPUTER) 


Case  I,  Satellite  User 


Number 
of  Navsat s 

Choice 

Number  of 
Time  Steps 

CPU  Time 
(sec) 

Amo  un  t 
of  Core 
(bytes) 

Costs 
(MUs ) 

24 

All 

48 

223 

300  K 

27.24 

24 

Zeni t h/nad i r 

48 

21 

300  K 

2.68 

3b 

All 

57 

348 

722  K 

42.70 

3 b 

Zenitli/nadir 

120 

b7 

300  K 

8.  31 

(;ase  II,  Single  Earth-Based  User 


24 

A1  1 

144 

25 

300  K 

3.91 

..  -._J 

Zen i t h 

1 44 

21 

300  K 

2.69 

Case  111,  Global  Net  of  Earth-Based  Users 


24 

Zenith  i 

1228 

300  K 

149.29 

‘'"All" 

refers  to  using  all  navsat s 

in  view,  4 

at  a time, 

, in  all 

calculations;  zenlth/nadir  refers  to  using  the  navsat  most  directly 
above  or  below  as  one  of  the  4 in  all  calculations. 

^'^Cost  in  Machine  Units — a direct  measure  of  computing  cost  used 
at  Rand,  which  accounts  for  core,  CPU  time,  and  input /output  resources 
used  for  a given  run;  currently,  a Machine  Unit  (Mil)  costs  about  $1.25. 


3T 
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IV.  PROGRAM  INPUT  AND  OUTPUT 


Table  5 (pages  32-36)  is  a listing  of  the  input  data  set  for  a 
Case  I run  (using  the  navsat  most  nearly  above  or  below  the  user  in 
all  calculations)  and  its  resulting  output.  This  run  treats  the  user 
satellite  as  a very  high  apogee  ballistic  missile.  Table  6 (pages  37- 
43)  is  a listing  of  the  input  data  set  for  a Case  II  run  (using  all 
navsats  in  view,  taken  four  at  a time)  and  its  resulting  output. 

Table  7 (pages  44-55)  is  a listing  of  the  input  data  set  for  a global 
distribution  of  earth-based  users  with  its  output. 


-33- 


Table  5 (Continued) 


ORBITAL  ELEMENTS 


m 

o 

ARGP 

RASC 

INC 

ANCM 

PER 

1 

0.0 

0.  0 

0.0 

63.00 

0.0 

12.00 

2 

0.0 

0.0 

0.0 

63.00 

45.00 

12.00 

3 

0.0 

0.  0 

0.0 

63.00 

90.00 

12.00 

4 

0.0 

0.0 

0.0 

63.00 

135.00 

12.00 

5 

0.0 

0.0 

0.0 

63. 00 

180.00 

12.00 

6 

0.0 

0.0 

0.0 

63.00 

225.00 

12.00 

7 

0.0 

0.0 

0.  0 

63.00 

270.00 

12.00 

(i 

0.0 

0.0 

0.0 

63.00 

315.00 

12.00 

9 

0.0 

0.  0 

120. 00 

63.00 

15.00 

12.00 

10 

0.0 

0.0 

120.00 

63.  00 

60.  00 

12.  00 

1 1 

0.0 

0.  0 

120.00 

63.00 

105.00 

12.00 

12 

c.o 

0.0 

120. 00 

63.00 

150.00 

12.00 

13 

0.0 

0.  0 

120.00 

63.00 

195.00 

12. 00 

14 

0.0 

0.0 

120. 00 

63.00 

240.00 

12.00 

15 

0.0 

0.0 

120.00 

63.00 

235.00 

12.00 

16 

0.0 

0.  0 

120.00 

63.00 

330.00 

12.00 

1 7 

0.0 

0.0 

240.00 

63.00 

30.00 

12.00 

1 d 

0.0 

0.  0 

240.00 

63.00 

75.00 

12.00 

19 

0.0 

0.0 

240.00 

63.00 

120.00 

12.00 

20 

0.0 

0.0 

240. 00 

63.00 

165.00 

12.00 

21 

0.0 

0.0 

240. 00 

63.00 

210.00 

12.00 

22 

0.0 

0.0 

240.00 

63.00 

255.00 

12.00 

23 

0.0 

0.0 

240.00 

63.00 

300.00 

12.00 

24 

0.0 

0.0 

240.00 

63.00 

345.00 

12.00 

USER  SATELLITE  ORBITAL  ELEMENTS 

25  0.82  11.24  47.50  144.21  69,  76  U.05 

TOTAL  TIMEIMINI  » 990 

time  INCREMENTIMINI  = 10 

BEAMWIDTH  ANGLEIOEGI  * 180.00 

FRACTICN  OF  NAVSAT  SPHERICAL  AREA  » 0.340 

THE  SATELLITE  MOST  NEARLY  ABOVE  OR  BELOW  IS  USED  AS  ONE  OF  THE  FOUR 
IN  ALL  CALCULATIONS  OF  THE  VOLUME  OF  THE  TETRAHEDRON 

^During  final  preparation  of  this  report,  the  authors  were  informed 
by  the  NAVSTAR/GPS  Joint  Program  Office  that  the  planned  orbital  con- 
figuration for  GPS  had  been  changed.  Users  who  wish  to  simulate  the 
currently  planned  configuration  should  interchange  the  anomalies  of 
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Satellite  user  calculations  terminate  when  and  if  the  orbit  intersects  the  earth 
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Table  6 (Continued) 


I 


I 


0KB1T4L  ELEMENTS 


ECC 

ARGP 

RASC 

INC 

ANOM 

PER 

1 

0.0 

0.  0 

0.0 

63.00 

0.0 

12.00 

2 

0.0 

0.0 

0.0 

63.00 

45.00 

12.00 

3 

0.0 

0.  0 

0.0 

63.00 

90.00 

12.00 

4 

0.0 

0.0 

0.0 

63.00 

135.00 

12.00 

5 

0.0 

0.  0 

0.  0 

63.00 

180.00 

12.00 

6 

0.0 

0.  0 

0.0 

63.00 

225.00 

12.00 

7 

0.0 

0.0 

0.0 

63.00 

270.00 

12.00 

U 

0.0 

0.0 

0.  0 

63.00 

315.00 

12.00 

y 

0.0 

0.0 

120.00 

63.00 

15.00 

12.00 

10 

0.0 

0.0 

120.00 

63.00 

60.00 

12.00 

1 1 

0.0 

0.0 

120. 00 

63.00 

105.00 

12.00 

12 

0.0 

0.0 

120.00 

6 3.00 

150.00 

12.00 

13 

0.0 

0.  0 

120. 00 

63.00 

195.00 

12.00 

14 

0.0 

0.0 

120. 00 

63.00 

240.00 

12.00 

15 

0.0 

0.  0 

120. 00 

63.00 

285.00 

12.00 

16 

0.0 

0.0 

120. 00 

63.00 

330.00 

12.00 

1 f 

0.0 

0.0 

240.00 

63.00 

30.00 

12.00 

18 

0.0 

0.0 

240. 00 

63.00 

75.00 

12.00 

19 

0.0 

0.0 

240. 00 

63.00 

120.00 

12.00 

20 

0.0 

0.0 

240.00 

63.00 

165.00 

12.00 

21 

0.  0 

0.  0 

240.00 

63.00 

210.00 

12.00 

22 

0.0 

0.  0 

240.00 

63.00 

255.00 

12.00 

23 

0.0 

0.0 

240. 00 

63.00 

300.00 

12.00 

24 

0.0 

0.0 

240.00 

63.00 

345. 00 

12.00 

USER  LOCATION  ON  EARTH 
LATITUDE  » U.O  DEGREES 
longitude  = 0.0  DEGREES 

MASKING  ANGLE  = 5.00  DEGREES 


I 

I 

I 


TOTAL  TIMEIMINI  = 
TIME  INCREMENTIMINI 


720 


ALL  SATELLITES.  TAKEN  FOUR  AT  A TIME.  ARE  USED  IN 
THE  CALCULATIONS  OF  THE  VOLUME  OF  THE  TETRAHEDRON 
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Table  7 (Continued) 
ORBITAL  ELEMENTS 


ECC 

ARGP 

RASC 

INC 

ANOM 

PER 

1 

0.0 

0.0 

0.0 

63.00 

0.0 

12.00 

2 

0.0 

0.0 

0.  0 

63.00 

45.  00 

12.00 

T 

0.0 

0.0 

0.0 

63.00 

90.00 

12.00 

4 

0.  0 

0.0 

0.0 

63.00 

135.00 

12.00 

5 

0.0 

0.0 

0.  0 

63.00 

IbO.OO 

12.00 

6 

0.0 

0.0 

0.0 

63.00 

225.00 

12.00 

7 

0.0 

0.  0 

0.0 

63.00 

270.00 

12.00 

U 

0.0 

0.0 

0.0 

63.00 

315.00 

12.00 

0.0 

0.  0 

120. 00 

63.00 

15.00 

12.00 

10 

0.0 

0.0 

120. 00 

63.  00 

60.00 

12.00 

1 1 

0.0 

0.0 

120.00 

63.00 

105.00 

12.00 

12 

o.c 

0.  0 

120. 00 

63.00 

150.00 

12.00 

13 

0.0 

0.  0 

120. 00 

63.  00 

195.00 

12.00 

14 

0.0 

0.  0 

120.00 

63.00 

240.00 

12.00 

15 

0.0 

0.  0 

120. 00 

63.00 

285.00 

12.00 

16 

0.0 

0.0 

120.00 

63.00 

330.00 

12.00 

1 7 

0.0 

0.  0 

240. 00 

63.00 

30.00 

12.00 

lb 

0.0 

0.0 

240.00 

63.  00 

75.00 

12.00 

19 

0.0 

0.0 

240.00 

63.00 

120.00 

12.00 

20 

0.0 

0.0 

240. 00 

63.00 

165.00 

12.00 

21 

0.0 

0.0 

240.00 

63.00 

210.00 

12,  00 

22 

0.  0 

0.  0 

240. 00 

63.00 

255.00 

12.00 

23 

0.0 

0.0 

240. 00 

63.  00 

300. 00 

12.  00 

24 

0.0 

0.0 

240.00 

63.00 

345.00 

12.00 

GLOBAL  DISTRIBUTION  CALCULATIONS 

MASKING  ANGLE  = 5.00  DEGREES 

LATITUDE  STEP  = 5.00  DEGREES 

LATITUDE  INCREMENT  = I 
LONGITUDE  INCREMENT  = I 

DILUTION  OF  PRECISION  PARAMETERS  PRINTED  AT  TIME  INCREMENT  OF  0 


TOTAL  TIMEIHINI  = 250 

TIME  INCREMENTIMINI  = 250 


1 


the  satellite  most  nearly  overhead  is  used  as  one  of  the  four 
IN  ALL  calculations  OF  THE  VOLUME  OF  THE  TETRAHEDRON 
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Appendix  A 

NAVSAT  SELECTION  ALGORITHMS  FOR  A SATELLITE  USER 


This  appendix  describes  the  algorithms  used  to  select,  from  all 
of  the  navsats  in  view  of  a satellite  user,  those  most  likely  to  yield 
small  DOP  values.  These  algorithms  were  developed  to  reduce  computa- 
tion costs,  accomplished  by  selecting  a value  for  the  input  parameter 
C (0  < C 1),  which  is  the  fraction  of  the  area  of  a sphere  with  radius 
equal  to  the  orbital  radius  of  the  navsats.  If  the  navsat  system  orbital 
configuration  involves  elliptic  orbits,  these  algorithms  cannot  be  used 

■k 

and  C should  be  set  equal  to  one.  Also,  computation  costs  can  be  re- 
duced further,  with  little  sacrifice  in  position  determination  accuracy, 
if  the  zenith/nadir  satellite  selection  mode  is  used  to  reduce  computa- 
tion time. 

If  the  navsat  orbital  configuration  involves  more  than  one  circular 
orbit  altitude  (a  "multishell"  system)  and  the  satellite  user  is  above 
all  of  the  shells,  the  algorithm  can  be  used  as  presented  below.  If  the 
satellite  user  is  between  or  below  the  shells  part  of  the  time,  the  situ- 
ation is  more  complex  because  there  are  two  possible  values  for  0^  (see 
below) — one  for  the  "zenith"  navsat  on  the  same  side  of  the  earth  as  the 
user,  and  one  for  the  "nadir"  navsat  on  the  other  side  of  the  earth 
from  the  user.  In  this  case,  suitable  values  of  C can  be  obtained  by 
careful  consideration  of  the  user/navsat  system  relative  geometries. 

As  with  elliptic  orbit  navsats,  C can  be  set  equal  to  one,  and  the 
zenith/nadir  navsat  selection  mode  can  be  used  to  minimize  computation 
time. 

The  following  is  a description  of  the  algorithms  for  use  with  any 
satellite  user  orbit,  as  long  as  the  navsat  system  uses  a single  altitude, 
circular  orbit.  The  number  of  navsats  considered  for  use  after  the 
zenith/nadir  navsat  is  selected  by  C;  the  larger  C is,  the  larger  the 
number  of  navsats  considered.  Values  of  C from  0.4  to  0.6  generally 
yield  good  results. 

*The  algorithms  can  be  used  if  the  satellite  user  is  always  above 
or  below  all  of  the  navsats  and  the  value  of  R,  which  Is  used  In  the 
expression  to  obtain  Oj,  Is  representative  of  the  navsat  orbits.  (For 
example,  set  R equal  to  the  average  of  the  semimajor  axes  of  the  navsat 
orbits.)  This  would  require  a minor  modification  of  the  computer  pro- 
gram. 


i 
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Tht>  two  algorithms  derived  below,  one  for  the  satellite  user  above 
the  navsats  and  the  other  for  the  satellite  user  below,  use  the  same 
criteria  for  selecting  the  first  navsat  of  the  four  to  be  used  for  the 
computation  of  the  DOPs , i.e.,  the  satellite,  of  those  in  view  of  the 
user,  which  is  closest  to  the  user's  zenith  or  nadir.  However,  the 
criteria  for  selecting  the  remaining  three  satellites  are  not  the  same 
for  the  two  algorithms. 


SELECTION  OF  THE  FIRST  NAVSAT 

For  convenience  of  explanation.  Fig.  A-1  shows  the  user  above  the 
navsats,  the  spherical  surface  containing  the  navsats,  and  the  earth 
plus  200  n mi  of  atmosphere  (hereafter  called  earth).  Also,  the  shaded 
area  that  is  shown  will  contain  the  navsats  which  are  retained  for  the 
computation  of  the  DOPs,  as  we  will  describe  later. 

At  each  time  point,  the  angles  6^  and  6^  are  computed  for  each 
user  and  navsat  combination  from 


-1  / R • P \* 


3 = cos 

n 


-1  / -II  • P 


If  0^  > if,  where  $ = cos  * ^ (f)’ 


then  the  satellite  is 


obscured  from  the  user's  view  by  the  earth  and  is  discarded  from  furtlier 
consideration  until  the  next  step. 

If  the  angle  a,  whicli  is  equal  to  n - (G  + B ) , is  greater  than 

n n 

the  navsat  antenna  beam  half-angle,  then  the  user  satellite  will  not  be 
illuminated  by  the  navsat  antenna  beam  and  the  navsat  will  not  be  used 
in  the  calculations. 

The  angle  Bj , corresponding  to  the  first  satellite  not  obscured  by 
the  earth,  and  the  number  of  the  satellite  are  stored.  Subsequent  B's 
are  computed  and  compared  to  the  B in  storage,  and  the  satellite  yield- 
ing the  smallest  B is  the  first  of  the  four  navsats  used  for  the  compu- 
tation of  DOP  values. 


P = |p|  (the  length  of  P)  and  R = |r|  (the  length  of  R) . 


satellite 


Fig.A-1 — Satellite  user  above  the  navsat  altitude 
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Kor  the  axis  of  the  navsat  antenna  pointing  toward  ttie  renter  of 
the  earth,  it  is  clear  that  any  navsat  above  the  plane  represented  by 
line  A-B  in  Fig.  A-1  must  have  an  antenna  beam  with  half-angle  greater 
than  90  deg  in  order  to  illuminate  the  satellite  user.  If  the  navsat 
is  below  the  plane,  the  antenna  beam  half-angle  required  will  be  less 
than  90  deg. 

SATELLITE  USER  ABOVE  THE  NAVSAT  SYSTEM 

If  the  remaining  three  navsats  needed  for  the  GDOP  computations 
are  on  or  near  the  tangent  circle  A-B  (Fig.  A-1)  and  are  separated  in 
longitude  by  120  deg,  then  the  volume  of  tlie  tetrahedron  formed,  as 
shown  in  Fig.  A-2  (a,h,c,d),  will  he  a near  maximum  and  the  correspond- 

■k 

ing  POOP  values  will  (almost  always)  be  minimum.  The  probability  that 
the  three  satellites  will  be  so  arranged  is  small;  thus,  all  of  the 
satellites  in  a band  (similar  to  a band  between  two  parallels  of  lati- 
tude) containing  the  tangent  circle,  and  whose  width  corresponds  to  the 
sum  of  the  central  angles  and  iS^,  are  retained  for  tlie  selection  of 
the  three  remaining  satellites.  The  satellites  not  in  this  band  are 
discarded  from  further  consideration  until  the  next  time  point.  The 
angles  and  are  obtained  from 

fS  = sin  ^ 0 > C , (A-2) 

1 \ sin  / T 

where  cos  0,^,  = |r|  / |p|,  and  ; or 

<5,  = 0,_  if  sin^  0„,  < C , (A-'i) 

and 

^2  = cos"^  - 2c|  - . 

The  central  angle  0,^  = cos  ^ j is  shown  in  Fig.  A-1.  The  quan- 
tity 0 < C < I (an  inptit  quantity)  represents  the  ratio  of  the  area  of 

* 

The  tangent  circle  A-B,  between  a cone  with  apex  at  the  user's 
position  and  the  spherical  surface  containing  the  navsats,  corresponds 
to  the  maximum  value  B can  have. 


'V'’ 
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tlu'  bund  to  thi“  area  of  the  spliere  containing  the  navsats.  If  the 
satellites  are  uniformly  distributed,  the  number  of  navsats  in  the  band 
will  be  the  product  of  C times  the  total  number  of  navsats.  Thus,  for 
C = 0.25,  the  number  of  CPS  satellites  in  the  band  (on  the  average) 
would  be  6.  Since  the  navsats  are  not  usually  uniformly  distributed, 

C S50.5  should  ensure  that  there  will  be  an  adequate  number  of  navsats 
available. 

SATinjJjrE_USjy<  BELOW  THE  NAVSAT  SYSTEM  ALTITUDE 

The  characteristics  of  the  algorithm  for  satellite  selection  for 
a satellite  user  below  the  navsats,  as  shown  in  Fig.  A-3,  are  not 
readily  apparent  if  it  is  desirable  to  reduce  the  computational  effort 
that  results  when  all  satellites  in  view  are  taken  four  at  a time.  Tlie 
number  of  different  combinations  of  four  satellites  can  be  large.  For 
example,  a user  near  the  altitude  of  CPS  will,  on  the  average,  see  22 
of  the  24  CPS  satellites,  and  the  number  of  different  combinations  of 
f our  will  be  7 315. 

The  satellite  selection  algorithm  presented  here  and  used  in  the 

program  is  based  on  the  simplifying  assumptions  that  the  first  of  the 

four  satellites  to  be  used  for  computing  CDOP  values  is  the  one  nearest 

the  user's  zenith  (or  nadir)  and  that  the  other  three  satellites  are 

in  or  near  a plane  which  is  normal  to  the  vector  P.  The  additional 

program  statements  required  to  select  the  optimum  orientation  of  the 

plane  were  not  used  because  of  singularity  problems  that  could  occur 

* 

when  the  user  is  just  below  the  navsat  altitude. 

Figure  A-3  shows  the  orbit  containing  the  first  of  the  four  nav- 
sats chosen  for  the  tetrahedron  volume  computation  and  the  band,  on 
the  spherical  surface  containing  the  navsats,  which  will  be  used  to 
select  the  other  three  navsats. 

The  angles  and  (?2  are  determined  as  before,  except  0.^.  is  now 
computed  from 

0 = p + cos  ' |— 1.^1—  cos  p|  , 

^ I |K|  I 

•k 

If  obtaining  minimum  CDOP  values  is  essential,  it  is  recommended 
that  the  algorithm  be  ignored  and  that  all  possible  combinations  of 
four  satellites  be  examined. 


i 


wliert'  p is  r lu*  angle  tn-tween  the  user's  horizuntal  and  a line  from  the 
usc-r  to  the  ein  le  repri'sen t ed  by  the  line  A-B  in  Fig.  A-3.  If  the 
first  navsat  chosen  is  on  the  user's  vertical,  tlicn  it  can  be  shown 
that  the  optimum  value  for  p is  about  19.5  deg.  This  value  of  p will 
maximize  the  volume  of  the  tetrahedron  formed  by  joining  the  tips  of 
the  four  unit  vectors  which  point  from  the  user  to  the  four  navsats. 

If  the  first  chosen  navsat  is  above  the  user's  horizontal  (as  shown  in 
the  figure),  the  angle  p will  be  below  the  user's  horizontal — and  vice 
versa . 

Figure  A- 1 shows  tliat  the  half-angle  of  the  navsat  antenna  beam, 
which  is  directed  toward  the  center  of  tlie  earth,  must  be  as  large  as 
the  angle  * if  the  user  satellite  is  to  receive  ranging  information 
from  a navsat  located  at  point  A. 
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Appendix  B 
PROGRAM  LISTING 


1. 

MAIN 

(pages  65-77 

2. 

ORBINI(I) 

[page  78] 

3. 

ORBIT(I,T,R,VEL,AC) 

[pages  79-80 

4. 

TRMATX(TR,I) 

[page  81 ] 

5. 

MATMUL(T,V,0) 

[page  82] 

6. 

POINT(ALO,ALA,TIM,VF.C) 

(page  83] 

7. 

TAAT(MAX,KXX, MATRIX) 

[page  84] 

8. 

VOLUME(OVEC, IDSAT,VOL) 

[page  85] 

9. 

COVNAV(G, ID,NAT,SIG) 

[page  86] 

10. 

BLOCK  DATA 

[page  87] 

LI. 

VECT0R(V1 ,I ,V2,V3) 

[page  88] 

12. 

D0T(V1,V2) 

[page  89] 

13. 

UNIVEC(V,UV) 

[page  90] 

14. 

TALL(KOT,KXX, MATRIX) 

[page  91] 
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HAIN 

COMMON/ORB IS /P 1 3 7, 25 ) /CON/C ( 1 8 ) 

DIMENSION  AA(3),G(37,A). lOS ( 37 ) , I S I C (4 ) . KXXI  58 905,4) ,NSGO(37) ,R(  3) 
1,RF(3) ,RHX(37, 3) . S I GT ( 6 ) ,U ( 3 ) , UPV( 3 ) ,USU V( 37 , 3) ,UTS ( 3) ,VE(3) ,Z(3), 
2 YEI3) .XNI  3) . ISAVEI  37) , I INT ( 3/ ) , UEL ( 37 ) , lOUT (37)  ,STN(3) 

DIMENSION  QS(<(36)  ,CAG(  19,18).GLEB(18),PI  B(  18),  ACL  ATI  19,36)  , 

1 ACT0T(36) ,0BLAT(19,36),0BDIS(36),CL(I9) ,MIN ( 19 , 36) , MA X( 19 , 36) , 

2 CAGXI 19,18) ,NSPL( 19) ,SKEG( 19 , 36,6 ) , SKEGXI 19 , 36,6 ) ,GLOS( 6 , 36 ) , 

3 COOP (19, 36, 6), I SAPP (32) ,PER(36) 

REAL  LAT , LONG, LATDEG,LADG, LAO 
KTR=48 

IPRINT=0 
Z( 1 )=0. 

Z(2)=0. 

Z(3)  = 1.E«-10 


LOC=l;  USER  ON  SATELLITE 

L0C=2;  USER  ON  GROUND  AT  SPECIFIED  LATITUDE  AND  LONGITUDE 
L0C=3;  GLOBAL  CALCULATIONS 


C 


C 


C 


READ  (5,175)  LOC  ,NJL  , I SCMP 
GO  TO  (100,110,120) , LOC 


100  DO  105  K*1 ,5 
READ  (5,180) 
105  CONTINUE 

READ  (5,180) 
NJM=NJL+1 
READ  (5,160) 
READ  (5,185) 
READ  (5,175) 
GC  TO  155 


(P(N,K) ,N^1,NJL) 

(PER(N) ,N=1,NJL) 

(P(NJM,K) ,K=1,5) 
PER(NJM) , AIN,CK 
INC, ITF 


110  DO  115  K«l,5 
READ  (5,180) 
115  CONTINUE 

READ  ( 5,180) 
READ  (5,185) 
READ  ( 5, 175) 
GO  TO  155 


(P(N,K) ,N=1,NJL ) 

(PER(N) ,N>1,NJL) 
ATL,ONGL,ELEVAT 
INC, ITF 


120  DO  125  K*1 ,5 

read  (5,160)  (P(N,0  ,N-1,NJL) 

125  CONTINUE 

READ  (5,180)  (PER(N) ,N-1,NJL) 

READ  (5,185)  LATOEG, ELEVAT 

READ  (5,175)  L AT IC , LONl C , INC , I TF , I PFREO, IT  I ME 

DO  135  MA»1,19 

CL(MA)«0. 

NSPL(MA)*0 
DO  130  MB>1,36 
MINIMA, MB) ^30 
MAXIMA, M6)>0 
DO  130  MC=l,6 
SKEGX(MA,MB,MC)=0. 

C0OP(MA,MB,MC)>0. 

130  SKEG(MA,MB,MC)*0. 

DO  135  MD= 1,  18 
CAGX(MA,MD )=0. 


o o o 
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MA  IN 

135  CAGIMA.MDI =0. 

DC  140  MA=1,36 
DO  140  MB=1,6 
140  GL0S(MB,MA)=0. 

DO  150  MA«l, 19 
ACT0T(MA)«0. 

OBDISIMAMO. 

DO  145  MD=1,36 
ACLAT(MA,MD)=0. 

145  0BLAT(MA,MD)=0. 

150  GLEB(MA)=0. 

C 

155  CONTINUE 

WRITE  (6.195) 

WRITE  (6,190) 

DO  160  IP=1,NJL 

WRITE  (6,215)  (IP,  (P(IP, IV) , IV  = l,5)  ,PER(  IP) ) 

160  CONTINUE 

IF(LOC.EO.l)  WRITE  (6,200) 

IP=NJL»l 

IF(LOC.EO.l)  WRITE  (6,215)  ( I P , ( P ( I P , I V) , I V= 1 , 5 ) , PER ( I P ) ) 

1FIL0C.E0.2)  WRITE  (6,210)  ATL ,ONGL ,ELEV AT 

IF(LOC.E0.3)  WRITE  (6,205)  EL  £ VAT , L ATDEG ,L AT IC , LON IC , I T I HE 

ITT=( ITF-1 )*INC 

WRITE  (6,220)  ITT, INC 

IF(LOC.EO.I)  write  (6,225)  AIN.CK 

IF(L0C.EQ. l.ANO. ISCMP.EQ.Ol  WRITE  (6,166) 

IF(LOC.NE. 1. AND. ISCHP.EQ.O ) WRITE  (6,165) 

IF(  ISCMP.EO.l ) WRITE  (6,170) 


THE  FOLLOWING  FORMATS  HAVE  TO  DO  WITH  INPUT 


165  FORMAT) IHO.lOX, 'THE  SATELLITE  MOST  NEARLY  OVERHEAD  IS  USED  AS  ONE 
lOF  THE  F0UR'/,11X, ' IN  ALL  CALCULATIONS  OF  THE  VOLUME  OF  THE  TETRAH 
2E0R0N* ) 

166  F0RMAT(1H0,10X,'THE  SATELLITE  MOST  NEARLY  ABOVE  OR  BELOW  IS  USED  A 
IS  ONE  OF  THE  F OUR* / , 1 1 X , • I N ALL  CALCULATIONS  OF  THE  VOLUME  OF  THE 
2TETRAHEDR0N' ) 

170  FORMATIIHO.IOX,' ALL  SATELLITES,  TAKEN  FOUR  AT  A TIME,  ARE  USED  IN' 
1,/,11X,'THE  CALCULATIONS  OF  THE  VOLUME  OF  THE  TETRAHEDRON') 

175  F0RMATI10I5) 

180  FORMAT! 12F6.0) 

185  FORMAT! 7F10.0) 

190  FORMAT! IHO, 2 2X,'  E CC • , 5X , ' APGP ' , 4X , ' RASC ' , 4X , ' I NC ' , 5X , 'ANOM' ,4X, 
I'PER' ,//) 

195  FCRMATdHl  ,//32X, 'ORBITAL  ELEMENTS') 

200  F0RNAT(IH0,10X,'USER  SATELLITE  ORBITAL  ELEMENTS'/) 

205  FORMATIIHO.IOX, 'GLOBAL  DISTRIBUTION  CALCULATIONS',//, 

1 IIX, 'MASKING  ANGLE  = ',F6.2,'  DEGREES',/, 

2 IIX, 'LATI  TUOE  STEP  ',F6.2,'  DEGREES',/, 

3 IIX, ' LATI TJDE  INCREMENT  = ',13,/, 

4 IIX, "LONGITUDE  INCREMENT  = ',13,/, 

5 IIX, ' DILUTION  OF  PRECISION  PARAMETERS  PRINTED  AT  TIME  INCREMENT  0 
6F  *,I5) 

210  FORMAT ( IHO, lOX, 'USER  LOCATION  ON  EARTH', /IIX, 

1 'LATITUDE  = *,F5.2,'  OEGRE E S ' / , I IX , ' LONG  I TUDE  * ',F5.2. 

2 • DEGREES',/,  IIX, 'MASKING  ANGLE  = ',F6.2,'  DEGREES') 

215  FORMAT(lH  , I 5X  , 1 3,  I X , i,F  8 . 2 ) 


non 


I 


1 


MAIN 

220  FORMATUHO.IOX,' TOTAL  TIHE(MIN)  = ',I5,/,ilX, 

I 'TIME  INCREMENT(MINJ  = ',13! 

225  FORMATdH  , I OX, ' BE  AMW  I OTH  ANGLE<DEG»  = • ,F6.  2 , / . I IX, 
I 'FRACTION  OF  MAVSAT  SPHERICAL  AREA  ',F6.3) 

SET  UP  ORBITAL  ELEMENTS 


( 


00  230  N*l ,NJL 

PIN,22)=22  808.*(  PERINI  Z2'*.  )♦♦(  2./3.  I 
CALL  ORBINI  (N) 

230  CONTINUE 

IF(L0C.E0.2.0R.L0C.E0.3»  GU  TO  235 
P( IP,22)=22b0d.* (PER (IP»/24.)**C2./3.» 

CALL  ORBINI  (IP) 

235  CONTINUE 
INCA=0 
NST0=0 
MAXNSS=0 
00  625  IT=1,  ITF 
TID=FLOAT( INCAI/14A0. 

INCA=INCA+1NC 
NN  = I P 

IF(LOC.E0.2.OR.L0C.E0.3I  NN=NJL 
DO  2'tO  N = 1,NN 

CALL  ORBIT  ( N, TI D , PER , RF , V E , A A ) 

DO  240  IV=l,3 
RMX(N,IV)=RF( IV) 

240  CONTINUE 

IF(L0C.E0.2I  GU  TO  335 
IF(L0C.E0.3)  GO  TO  330 
C 
C 

DO  245  IV=l,3 
UPV(IV)=RMX(  IP, IV) 

245  CONTINUE 
IK=0 

BETAMN=«4. 000000 
DO  265  N»1 ,NJL 
DO  250  IV*l,3 
R( IV)=RMX( N, IV) 

250  CONTINUE 

CALL  VECTOR  ( R,2 ,UP V , UTS ) 

AR*SQRT( 00T( R,R) ) 

AP»SQRT(OOT(UPV,UPV) ) 

A^ARCOS ( 2. 0926144E*07/AR ) 

IF(AP.LE.2.0926I44E«^07)  WRITE  (6,251 ) 

251  FORMATdH  , / ,3X,  • TERM  I NA  T I ON  OF  RUN,  ALTITUDE  APPROCHING  ZERO'! 
IF(AP.LE.2.0926144E«^07)  STOP 
B=ARCUS(2.0926144£»07/AP) 

PHdA^B 

OOTRP=iOOT(  R,UPVI 
THETN«ARCOS(OOTRP/ (AR*AP) ) 

IF(THETN.GE.PHI)  GO  TO  265 
IK-IK»l 
ISAVEI IK)>N 
AU-SORT(OOT(UTS,UTSI ) 

OOTUP»OOT(UTS,UPV) 

BETAN«ARCOS(-OOTUP/ ( AU*AP)  ) 


J 

1 
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MMN 


AIWID=A1N/C(2) 

BWIDTH=C  (3)-(THETN4-BETAN) 

IF( AIWID.LT.BWIDTH)  go  to  265 
BETASV=BETAN 

IF<BETAN.GT.C(5I » BETAN=CC 3)-BETAN 
255  IF<BETAN.LT.BETAMN»  GO  TO  260 
GO  TO  265 
260  ISATNO^U 

BETAMN*BETAN 
BETAS=BETASV 
265  CONTINUE 

DO  285  I*l  ,IK 

IF(  ISAVEU  ) .EO.ISATNO)  GO  TO  270 
GO  TO  285 
270  INDX«I 

DO  275  KI=1, INDX 
KII=KIH 

IINTIKI I 1= ISAVFIKI ) 

275  CONTINUE 

IN0X=1N0X>-1 
DO  280  KI-INDX.IK 
IINTIKI »=ISAVE(KI» 

280  CONTINUE 
GO  TO  290 
285  CONTINUE 
290  IINTm  = ISATNO 

IF! ISCMP.EO. n GO  TO  335 

THE  NUMBERS  OF  THE  SATELLITES  WHICH  FIT  CRITERIA  FOR  USE,  PLUS  THE 
ONE  NEAREST  TO  OVERHEAD  HAVE  BEEN  CALCULATED 

lOSI 11  = 1 INTI  II 
NSS=l 

DO  295  IV=1,3 
RUVI=RMX(  IINTII),1VI 
295  CONTINUE 

CALL  VECTOR  I R, 2 ,UP V , UTS  I 
CALL  UNIVEC  lUTS.U) 

DO  300  IV=l,3 

USUVI I INTI  1 1 ,IV)=UI 1 V) 

300  CONTINUE 

DO  325  N=2,IK 
DO  305  1V=1,3 
RIIVI=RMX( I INTINI, IVI 
305  CONTINUE 

CALL  VECTOR  IR ,2 ,UPV , UTS  I 
CALL  UNIVEC  IUTS,U) 

DO  310  IV=1,3 
USUVIIINTIN)  ,IV)=UIIVI 
310  CONTINUE 

AR*SQRTIDOTIR,RH 
AP=SORTIOOTIUPV,UPVI I 
IFIAP.GT.ARI  THETT=ARCJSIAR/AP) 

IFIAP.lE.AR.AND. BETAS.GI .C I 5 I I THET T= 1 09 . 5/C  I 2 l-ARS I NI AP* . 33 3806 86 
l/ARI 

IFI AP.LE.AP.ANO.BETAS.LE  .Cl  51  I THFTT-70. 5/C  I 2 l-ARS IN  I AP*. 3 3380686/ 
lAR) 

DOTRP=OOTI R.UPVI 


— "V 
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THETN»ARC0S(D0TRP/ (AR*API ) 

SINTT=SIN( THETTI 
SlNTSa=SINTT*'»2 
IFISINTSQ.GT.CKI  GO  TO  315 
OELONE=THETT 
ANG*l.-(CK/.5» 

DELTW0=*ARC0S(ANGI-THETT 
GO  TO  320 
315  CONTINUE 

OELONE=ARS INICK/SINTT  ) 

DELTWO=OELONE 

320  IFITHETN.LT. (THETT-DELONEI » GO  TO  325 
IF(ThETN.GT.(THETT^OELTWJ> J GO  TO  325 
NSS  = NSS+  I 
IOS(NSS»=I INT(N) 

325  CONTINUE 

IF(NSS.LE.4)  GO  TO  625 
GO  TO  335 

END  OF  USER  ON  SATELLITE 

CALCULATIONS  FOR  USER  ON  EARTH  AT  SPECIFIED  LAT  AND  LONGt  AND  GLOBAL 
DISTRIBUTIONS  FOLLOW 
330  CONTINUE 
LKL=36 
LLL*19 

335  CONTINUE 

IFILOC.EQ.  2.0R.L0C.EQ.H  LKL=1 
IF(L0C.E0.2.0R.L0C.EQ.l)  LLL»l 
DO  500  K=l ,LKL ,LONIC 
DO  500  L=1  ,LLL,LATi: 

IFILOC.EQ. l.ANO.ISCHP.EO.O)  GO  TO  <.00 
IFIlOC.EO. 1. AND. ISCHP.EO.l 1 GO  TO  375 
LONG=FLOAT (K-1 JAIO. 

LAT*90.-FL0AT(L-n  *LATDEG 
IFILOC.EO. 2)  LONG=OnGL 
IFILOC.EO. 2)  LAT=ATL 
CALL  POINT  ILONG.LAT ,T ID.UPVI 
NSS=0 

CHL»»CLIL)H. 

DO  350  N>1 «NJL 

00  340  IV- I, 3 
RI IVI-RMXI N. IV) 

340  CONTINUE 

CALL  VECTOR  lUPV ,2 t R , STN I 

SE— OOTISTN,UPV)/SORTIDOTISTN,STN)*DCTIUPV,UPV) ) 

1 FI ABSI SE) .GE.. 9999999)  SE=S IGNI 1.  ,SE ) 

EL-ARSINISE)*CI2) 

IFIEL.lt. ELEVAT)  GO  TO  350 
NSS-NSS+l 
I INTINSS )-N 
UELIN)-EL 

CALL  VECTOR  I R , 2 .UP V , UT S ) 

CALL  UNIVEC  lUTS.U) 

00  345  IV-l,3 
USUVIN,IV)-UIIV) 

345  CONTINUE 
350  CONTINUE 


ooooonoo 
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IF(NSS.GE.A)  GO  TO  355 
IFINSS.lt. 4»  GO  TO  500 
355  CONTINUE 

NSTO=NSTO«-NSS 
NSPL(L)»NS  PL  (U^NSS 
NSP=NSS*1 

MAXNSS=>MAXO(  HAXNSS.NSSi 
HIGH=0. 

00  365  NUM=1,NSS 
IX=1INT(NUM) 

AEA  = UEL I IX  ) 

IF(AEA.GT.HIGH)  go  to  360 
IFIAEA.L E.HIGHI  GU  TO  365 
360  HIGH*AEA 
NN^NUM 

365  CONTINUE 

DO  370  NU= 1 NSS 
IFINU.EQ.NNI  NX-l 
IF(NU.LT.NN)  NX=NU*1 
IFINU.GT.NNI  NX*NU 
iaS(NX»«IINT(NU» 

370  CONTINUE 
GO  TO  ^00 

lOSINSS)  HAS  the  SATELLITE  NEAREST  TO  OVERHEAD,  THEN  ALL  OTHERS 
WHICH  FIT  CRITERIA 

NEXT  CALCULATIONS  ARE  CONCERNED  WITH  FINDING  THE  COMBINATION  OF 
FOUR  SATELLITES  WHICH  HAVF  THE  GREATEST  VALUE  OF  THE  VOLUME  OF  THE 
TETRAHEDRON  FORMED  BY  THEM 

375  CONTINUE 

DO  390  IKN=l,IK 
DO  380  IV=l,J 
RnV)=RMXI  IINT(IKN).  IV) 

380  CONTINUE 

CALL  VECTOR  ( R , 2 , UPV , UTS  I 
CALL  UNIVEC  (UTS.U) 

DO  385  IV* 1,3 

USUVI I INTI  IKN) ,IV)=J{ IV) 

385  CONTINUE 
390  CONTINUE 

00  395  I JK=l , IK 
395  lOSI  1 JK)*I  INTI  UK) 

-tOO  B0X«-10. 

IFI ISCMP.E  0. 1)  GO  TO  405 

KOT«NSS 

NSS*NSS-l 

CALL  TAAT  I NSS , KCOM , KXX ) 

I SICI I )= lUSI  I ) 

LPN*2 
GO  TO  410 
405  KOT=IK 

IF  ILOC.E0.2.0R.LOC.E0.3)  KOT*NSS 
CALL  TALL  IKOT , KCOM, KXX ) 

LPN=I 

410  DU  430  M=1 ,KCOM 


) 

I 

j 

i 


1 

( 

i 


) 

I 


■ i 
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1 

( 
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MAIN 

DO  ^15  LPQ=LPN,'» 

NUX»KXX( M, LPU) 

ISICaPO)=  lOS(NUX) 

<il5  CONTINUE 

CALL  VOLUME  (USUV, I SI C , VOLUMJ 
IF(VOLUM.GT.BOX)  GO  TO  ^20 
IF(VOLUM.LE.BOX)  GO  TO  <f30 
<•20  BOX=VOLUM 

DO  MC=Vt<» 

NSGD(MC) =I SICIMCI 
425  CONTINUE 
430  CONTINUE 

DO  440  N = 1 ,4 
DO  435  IV=lf  3 
R(IV)=RMX( NSGDINI, IV) 

435  CONTINUE 

CALL  VECTOR  I R , 2 ,UP V , UTS  ) 

CALL  VECTOR  ( Z,3  ,UPV, YE) 

CALL  VECTOR  (UPV , 3 • YE , XN ) 
AX=SORT(DOT(XN,XN) ) 

AY=SORT(DOT( YE, YE)  ) 
AP=SURTIOOT(UPV,UPV) ) 

ALT  = AP/607  6. 1 16-3444 . 

AS=SORT(OOT( UTS.UTS) ) 

G(NSGDIN) , I )=00T ( XN, UTS ) / ( AX*AS ) 
GINSGDIN) ,2) =00: (YE, UTS) /(AY* AS) 
G(NSGD(N), 3)=D0T (UPV, JTS)/(AP*AS) 
G(NSGD(N) ,4) =l  . 

440  CONTINUE 

CALL  COVNAV  ( G ,NSGD , 4 , S I GT  ) 
IF(L0C.E0.3)  GO  TO  480 
IF(LOC.E0.2)  ALT=0. 

DO  445  11=1,4 
IOUT( I I )=NSGD(  II  ) 

445  CONTINUE 
I L = 1 
LK  = 4 

DO  455  NN=1,K0T 


I TST  = NSGO(  IL  ) 

IF( ITST.EO.IQS(NN) ) GO  TO  450 
LK=LK*1 

IOUT(LK)«=IOS(NN) 

GO  TO  455 
450  IL=IL*1 

IF(IL.EQ.5)  GO  TO  460 
455  CONTINUE 
460  K0S«NN*1 

DO  465  LK=KOS,KOT 


IOUT(LK)=I JS(LK) 

465  CONTINUE 

ITOUT*( IT-1)*INC 
KTR=KTR*2 

IFIKOT.GT.  18)  KTR=KTR*1 
IF(KTR.GE.50)  GO  TO  470 
GO  TO  475 
470  WRITE  (6,525) 

KTR  = 0 

475  WRITE  (6,530)  I T OUT  , AL T , ( S I GT  ( K P ) , KP= I ,6 ) , ( I OUT ( K P)  , KP»  1 , KOT ) 


o o o o r>  o 


END  Of  CALCULATION  FOR  SINGLE  USER  ON  EARTH 
GO  TO  625 

FOLLOWING  CALCULATIONS  FOR  GLOBAL  DISTRIBUTION 

400  DO  405  I00P=l,6 
485  CDOP(LfK,IOOP)=SIGT( IDOP) 

DO  490  NS-1,6 

KA:*MAXO(  1,MIN0UNT(SIGTCNSI*5.»1.I  ,361) 

490  SKEG(LfKA,NS)=SKEG(LtKA,NS)«-l. 

DO  495  NA=l,KOT 
I-IOSINA) 

ELT*UEL(  I) 

KQ=MIN0(INTIELT/5.)«-l,lb) 

495  CAG(L,KO)=CAG( L,KOi*l. 

MAX(L.K) =MAX0(MAX(L,KI ,KOT» 

MIN(L,KJ«:MINO(MIN(l,K»  ,K0T) 

OBLAT(LtNSP) =OBLAT(LtNSPI>l. 

500  CONTINUE 

IFILOC.EO.B.AND.IPFREO.EO.O)  GO  TO  625 
C NO  INTERMEDIATE  PRINT 

IFIIPRINT.EQ.O.AND.IT.EO.n  GO  TO  505 
C PRINT  FIRST  TIME  STEP 

IPRINT»I PR  INT  + l 
IFIIPRINT. EO.ITIME)  GO  TO  505 

C PRINT  EACH  TIME  STEP  REQUESTED 

GO  TO  625 
505  IPRINT»0 


DO  515  I00P=1,6 

IFUOOP.EQ.n  WRITE  (6,535)  ITIME 
;F( ID0P.E0.2 ) WRITE  (6,540)  ITIME 
IF( I00P.E0.3)  WRITE  (6,545)  ITIME 
IF( ID0P.EQ.4)  WRITE  (6,550)  ITIME 
IF( I00P.EQ.5)  WRITE  (6,555)  ITIME 
IFIID0P.E0.6)  WRITE  (6,560)  ITIME 
IFILATDEG.EO.IO. ) WRITE  (6,565) 

IF(LATO£G.E0.5.)  WRITE  (6,570) 

ICT  — 10 

DO  510  IC=1,36 
ICT»ICT+10 

IF(ICT.LE.90.UR.ICT.GE.190)  WRITE  (6,575)  I CT , (CDOP( IK , I C , IDOP ) , IK 
1-1,19) 

IF( ICT.EO.lOO)  WRITE  (6,580)  ( COOP ( I K , IC , IDOP ) , IK-1 , 19 ) 

IF( ICT.EO.llO)  WRITE  (6,585)  ( COOP( IK , IC , IDOP ) , I K-1 , 19 ) 

IF( ICT.E0.120)  WRITE  (6,590)  ( CDOP ( IK , IC , IDOP  ) , I K-1 , 19 ) 

IF( ICT.EQ.130)  WRITE  (6,595)  ( CDOP( IK , IC , IDOP ) , I K-1 , 19  I 
IF(ICT.E0.140)  WRITE  (6,600)  ( CDOP ( IK , IC , I OOP ) , I K-1 , 19 ) 

IF( ICT.EO. 150)  WRITE  (6,605)  ( CDOP ( IK , IC , IDOP ) , I K-I , 19) 

IF( ICT.E0.160)  WRITE  (6,610)  ( CDOP ( 1 K , IC , 1 OOP ) , I K -1 , 19 ) 

IF( ICT.EO. 1 70)  WRITE 
IF(ICT.EO.iaO)  WRITE 
510  CONTINUE 
515  CONTINUE 

DO  520  ICL-1,19 
DO  520  ICK=I,36 
DO  520  ICD-1,6 
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HAIN 

520  COOP(1CL,ICK,ICD»=0. 

THE  FOLLOWING  FORMATS  HAVE  TO  DO  WITH  A SPECIFIED  TIME  STEP  REQUEST 
FOR  PRINTING 


I 

r 

i 


I 

I 

i 


I 

I 


c 


•,I6,20X,'TD0P 
• . I6,20X, ’POOP 


52  5 FORMAT (1  HI , / , I X , ' T I ME ( MN ) • , 1 X , • ALT ( NM ) • , 

1 5X, • HDUP*  ,5X, 'MDOP*  , 5X, ‘TDOP’ t 5X, • PDOP* 

2 'SATELLITES  CHOSEN',//! 

530  FORMAT(1HO,I7,F9.0,6(F9.3»,IX,  18(IX,I2), 
535  FORMAT ( IHl  ,//, lOX, ' TI ME  = ' , I 6 , 2 OX , ' VD OP 
540  FORMATIlHl ,// ,10X, ' T IME 
lONTAL  PLANE' ) 

545  FORMAT(lHi,//,lOX, 'TIME 
1 ITION  ERROR' ) 

550  FORMAT! IHl ,//,10X, 'TIME 
555  FORMAT! IHl ,//,lOX, ' TI ME 
IION  ERROR' ) 

560  FORMAT ! IHl ,// , lOX, ' T I ME 
ION  ERROR' 1 

565  FORMAT !1H0,52X,' LAT ITUDE ',// ,9X, ' 90', 3X 

1 3X,'  50', 3X,'  40', 3X,'  30', 3X,'  20', 3X 

2 3X, '-20' ,3X, '-30' , 3X, '-40' ,3X, '-50' ,3X 

3 3X,'-90',/I 

570  F0RMAT(IH0,52X, 'LATITUDE', //,9X,'  90', 3X 

1 3X,'  70', 3X,'  65', 3X,'  60', 3X,'  55', 3X 

2 3X,'  35', 3X,'  30', 3X,'  25', 3X,'  20', 3X 

3 3X,'  0',/) 

575  FORMATIIH  , 2 X , I 3 , 1 X , 19F 6 . 2 I 

580  FORMATIIH  , ' L'  ,IX, ' 100' , IX, 19F6.2) 

585  FORMATIIH  , ' O' , I X, ' U O' , IX , 1 9F 6. 2 I 
590  FORMATUH  , ' N' , I X,  ' I 20' , IX,  1 9F  6.  2 ) 

595  FC'MATIIH  , ' G ' , IX , ' 1 30 ' , IX , 1 9 F 6 . 2 I 
600  FORMATIIH  , ' I ' , 1 X , ' I 4 O' , 1 X , 1 9F 6. 2 ) 

605  FORMATIIH  , ' T ' , I X , ' 1 50* , I X , I 9F 6.2 ) 

610  FORMATIIH  , ' U' , IX, ' 1 60 ' , 1 X , 1 9F 6. 2 I 
615  FORMATIIH  , ' D' , 1 X , ' I 70 ' , IX , I 9F 6. 2 ) 

620  FORMATIIH  , ' E ' ,l X, ' I 80 ' , IX , 19F 6 . 2 ) 


4X, ' VDOP' , 

,5X, 'GOOP' ,4X, 

/.72X, I8I IX,  I2J  i 
- ALTITUDE* ) 

■ , I6.20X, 'HDOP  - POSITION  ERROR  IN 

' , I6,20X, 'MDOP  - LARGER  COMPONENT 


- TIME' ) 

- three  dimensional 


HORIZ 
CF  POS 

POSIT 


', 16.20X, 'GDOP  - FOUR  DIMENSIONAL  POS'TI 


3X,' 

70' 

,3X,' 

60' 

3X,' 

O' 

, 3X,  '■ 

-10' 

3X, '-70' 

,3X,'- 

-80' 

3X, ' 

80' 

,3X,  ' 

75' 

3X,' 

45' 

,3X,* 

40* 

3X,' 

10' 

,3X,  ' 

5' 

625  CONTINUE 

IFILOC.EO. 1.0R.L0C.EQ.2)  GO  TO  920 
MNSSP0=MAXNSS^1 
DO  630  L 1 = 1 ,36 
630  QSRILI )=FL0ATlLI-l)*.2 
DO  640  I X= 1, Id 
PIBI 1X1=0. 

DO  635  I V=l,  19 

635  PIB(IX)*PI8I  IXI^CAGI 1Y,IX) 

640  PIBI I X)= IP IBI IX)/FL3ATINST0) ) *100. 

DO  660  LA=l,  19 

IF(CLILA).EO.O..OR.NSPL(LA).EO.O»  GO  TO  660 
DO  650  N=1 ,36 
DO  650  1=1,6 
SKEGXILA,N,I)=0. 

DO  645  J=N,3b 

645  SKEGXILA  ,N,II=SKEGXILA,N,I)»SKEGILA,J,n 
65  0 SKEGX(LA,N  , I ) =SKESX I LA , N , I ) /CL IL A ) 

00  660  LC=l,  18 
CAGXIlA,LC )=0. 
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MAIN 


DO  65&  ; F=LC.  la 

655  CAGXILA.LC ) »CAGX (L A, LC ) ♦CAGI L A , LF  ) 

CAGX(LA,LC  »=CAGXlLA,LCl/FLOAT(NSPL  (LAU 
660  CONTINUE 

WRITE  (6,780) 

WRITE  (6,785)  ( P I B ( I X ) , I X= 1 , 1 8 ) 

IF(LAT0EG.EQ.5.)  WRITE  (6,795) 

IF(LATDEG.E0.10. ) WRITE  (6,800) 

ICT*-5 

DO  665  LC=l,l8 
ICT=ICT*5 

WRITE  (6,815)  I:T, (CAGX(LA,LC) ,LA=1, 19) 

665  CONTINUE 
CVS=0. 

DO  680  LM= I, 19 

IF(CL( LM). EQ.0..0R.NSPL (LM) .EQ.O)  GO  TO  680 
CV=COS( (90.-FIQAT(LM-1)*LATDEG)*C( 1) ) 

CVS=CVS*CV 
DO  670  1 N= 1, 6 
00  670  11=1,36 

67  0 GLOS( IN, II ) = GL0S( IN, I I H-SKEGX(LM, I I,  IN)*CV 
00  675  LN=1,  18 

675  GLEB(LN)=GLEB(LN)>CAGX(LM,LN)*CV 
680  CONTINUE 

00  685  LN=  1.  18 
685  GLEB(LN)=GLEB(LN)/CVS 
WRITE  (6,805) 

WRITE  (6,310)  (GLEB(iC),IC=l,18) 

DO  690  IN=1,6 
00  690  I 1*1, 36 

690  GLOS( IN, II l = GLOS ( IN, I I )/CVS 
00  695  JD0P=1,6 
IF( JOOP.EQ.l  ) WRITE  (6,830) 

IF( J00P.EQ.2 ) WRITt  (6,835) 

IF( JD0P.E0.3)  WRITE  (6,8^0) 

IF( J00P.E0.4)  WRITE  (6,845) 

IF( JD0P.EQ.5)  WRITE  (6,8501 
IF( J00P.E0.6)  WRITE  (6,855) 

IF(LATDEG.EQ.10.)  WRITE  (6,820) 

IF(LAT0EG.EQ.5.)  WRITE  (6,825) 

00  695  I0SR=^1,36 

WRITE  (6,775)  QSR( I QSR ) , (SKEGX ( L K , I QSR, J OOP ) , LK=  1 , 19 ) 
695  CONTINUE 

WRITE  (6,860) 

00  700  I(3SR  = 1,36 

WRITE  (6,865)  QSR ( I QSR ) , ( GLO S ( I N , I QSR ) , I N= 1 , 6 ) 

700  CONTINUE 

WRITE  (6,870) 

LA0G=90. 

IFILATOEG.EO.IO. ) LAC=-10. 

''■•(LAT0EG.EO.5.)  LAD=-5. 

WRITE  (6.880) 

00  705  IL=1,19 

WRITE  (6,885)  L A OG , ( M A X ( I L , 1 K ) , I K = 1 , 36 ) 

LADG=LADG»LAD 
705  CONTINUE 

WRITE  (6,875) 

L A0G=90. 
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MAIN 

OU  715  IL=1,  19 
DO  710  IX=l,36 

IF(MIN(IL,IX).E0.30»  MIN(IL,IX»=0 
710  CONTINUE 

WRITE  (6,8851  L ADG, ( M I N( I L , I K ) , I K= 1 , 36  I 
LAOG»LADG^LAO 
715  CONTINUE 

DO  720  N=1 .MNSSPO 
0BDIS(N)=0. 

DO  720  L*1 ,19 

IF(CL(U  .E0.0..0R.NSPL(L).E0.0I  GO  TO  720 
OBLAT(L,N) *(OBLAT(L,N»/CL(U )*100. 

720  CONTINUE 

00  730  L=l,19  \ 

DO  730  N=l, MNSSPO  | 

ACLAT(L,N»=0.  j 

DO  725  M=N, MNSSPO  j 

725  ACLAT(L,N)=ACLAT(L,N»«-OBLAT(L,M)  ( 

730  CONTINUE  'j 

co=o. 

DO  740  L=l ,19  1 

1 F«CL( U .EO.O..OR.NSPL (L ) .EQ.O  I GO  TO  740 
CA=COS(  ( 90. -FLOAT!  L-1  J^LATDEGI  *Cm  ) 

CO=CO+CA 

DO  735  N=l, MNSSPO  j 

735  OBDIS(N)*:OBDIS(N)«-OBLAT(L,N»*CA 
740  CONTINUE 

DO  745  N=1 , MNSSPO 
745  OOOIS(NI=OBDIS(N)/CO 
DO  750  N=1 , MNSSPO 
ACT0T(N)*0. 

DO  750  M=N, MNSSPO  1 

150  ACTOT(N)=ACTUT(N)«-OBDIS(M) 

DO  755  I =l  ,32 
I SAPPI I ) *1-1 
755  CONTINUE 
NP*1 

MNS=MNSSPO 

IF  (MNSSPO. GT. 16)  MNS*16 
WRITE  (6,890)  ( I SAPP ( I ) , I = 1 , 1 6 ) 

760  WRITE  (6,900) 

L ADG=90. 

DO  765  IL=l,19 

WRITE  (6,895)  LAD3,  (OBLAT( IL ,N) ,N=NP,MNS ) 

LADG=LADG*LAO 
765  CONTINUE 

WRITE  (6,905) 

LA0G=90. 

DO  770  IL* 1,  19 

WRITE  (6,895)  LAL/G , ( ACL  AT  ( IL  , N ) , N=  NP,  MNS  ) 

LADG*LA0G*LA0 

770  CONTINUE  i 

WRITE  (6,910)  . 

WRITE  (6,790)  ( OBD I S ( N ) , N = NP , MNS  ) j 

WRITE  (6,915)  I; 

WRITE  (6,790)  (ACT0T(N),N*NP,MNS)  I 

IFIMNSSPO.lt. 16)  GO  TO  920  i 

IFIMNS.GT.16)  GJ  TO  920  I 
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MNS=MNSSPG 

NP=17 

WRITE  (6,890»  (I SA PP < I) , I = 17 , 3 2 J 
GO  TO  760 

THE  FOLLOWING  FORMATS  HAVE  TO  DO  WITH  A GLOBAL  SYSTEM 


7 75  FORMATIIH  , F 4 . 1 , 3X , 1 9F6 . 3 ) 

780  FORMAT ( IHl , 7/7/7 IX, • ELEVATION  DISTRIBUTION  - PROBABILITY  THAT  THE 
ISATELLITE  IN  VIEW  WILL  HAVE  ELEVATION  ANGLES  AS  L ISTED*  , 7/ ,<f8X , 

2  'ELEVATION  ANGLE • ) 

785  FORMAT! IHO ,6X, • 0-5  5-10  10-15  15-20  20-25  25-30  30-35  35-40  40 

1-45  45-50  50-55  55-60  60-65  65-70  70-75  75-80  80-85  85-90* .//, 2X , 

2  'PROS' , IX, 18F6.1) 

790  FORMAT! IhO , 3X , ' PRQB ' , 5X , 16F7 . 2 » 

795  FORMAT! IHO, 77/7, IX, ' LATITUDE  ELEVATION  DISTRIBUTION',/, 

1 IX, • PROBABIL ITY  THAT  ANY  SATELLITE  IN  VIEW  WILL  HAVE  AN  ELEVAT 
2I0N  ANGLE  GREATER  THAN  OR  EQUAL  TO  THOSE  L ISTEO'  ,//, 52X, 

3't ATITUDE' ,7,2X,'ELEV',7.1X,'ANGLE*,  U , ' 90  85  80  75 

4 70  65  60  55  50  45  40  35  30  25  20 

5 15  10  5 0',7» 

800  FORMAT ! IHO ,77/7, IX, ' LATITUDE  ELEVATION  DISTRIBUTION',/, 

1 IX, ' PROBABIL ITY  THAT  ANY  SATELLITE  IN  VIEW  WILL  HAVE  AN  ELEVAT 

2I0N  ANGLE  GREATER  THAN  OR  EQUAL  TO  THOSE  L I S TED' , 7/ , 52X , 
3'LATITUDE' ,/,2X,'ELEV',/,lX,'ANGLE' , IX, • 90  80  70  60 

4 50  40  30  20  10  0 -10  -20  -30  -40  -50 

5 - 60  -70  -80  -90', 7) 

805  FORMATIIHO, 7///, IX,* ACCUMULATIVE  ELEVATION  DISTRIBUTION',/, 

IIX, 'PROBABILITY  THAT  THE  ELEVATION  ANGLE  IS  GREATER  THAN  OR  EQUAL 
2T0  THOSE  L ISTED' ,77, 48X, 'ELEVATION  ANGLE') 

810  FORMATIIHO, 7X,'  0 5 10  15  20  25 

10  45  50  55  60  65  70  75  80 

2 'PROS', IX, 18F6.1) 

815  FORMATIIH  , I 3,4X, 1 9F6.2 ) 

820  F ORMAT ! 1 HO, 52 X,' LATITUDE* ,/7,2X, *NUM',5X, 

1 ' 90', 3X,'  80', 3X,'  70', 3X,'  60', 

2 3X,'  50', 3X,'  40', 3X,'  30', 3X,'  20',3X,»  10', 3X,'  0',3X,'-10', 

3 3X, '-20' ,3X, '-30* , 3X,'-40* ,3X, '-50* ,3X,'-60* ,3X , '-70* , 3X, '-80' , 

4 3X,*-90',7) 

825  F0RMATI1H0,52X,' LAT ITUDE* ,77 ,2X, 'NUM* ,5X, 

1 


30  35 

85* ,77, 2X, 


2 3X,  ' 70*  , 3X, • 65 

3 3X, • 35*  ,3X,  ' 30 

4 3X,  ' O'  ,7  I 

830  FORMATIlHl ,777, IX, 
IISTRIBUTION* ,7, IX, 
2HAN  number  L ISTED* 
835  FORMAT! IHl ,777 , IX, 
IISTRIBUTION* ,7, IX, 
2NE  WILL  BE  GREATER 
840  FORMATIlHl ,777, IX, 
1 ISTRIBUTION* ,7, IX, 
20P  WILL  BE  GREATER 
845  FORMATIlHl ,777, IX, 
IISTRIBUTION*  ,7, IX, 
/NUMBER  LISTED',//) 
850  FORMAT! IHl ,7/7 ,1X, 
IISTRIBUTION* ,7, IX, 


• 90', 3X,'  85', 3X,'  80', 3X,*  75', 
,3X,'  60', 3X,*  55', 3X,'  50', 3X,'  45', 3X,*  40', 
,3X,'  25', 3X,'  20',3X,'  15»,3X,'  10‘,3X,*  5», 

DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  D 
PROBABILITY  THAT  ALTITUDE  DOP  WILL  BE  GREATER  T 
77) 

DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  D 
probability  that  POSITION  OOP  IN  HORIZONTAL  PLA 
THAN  number  LISTED', 77) 

DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  0 
probability  that  larger  COMPONENT  OF  POSITION  D 
THAN  number  LISTED', 77) 

DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  D 
PROBABILITY  THAT  TIME  OOP  WILL  BE  GREATER  THAN 

DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  0 
probability  that  three  DIMENSIONAL  POSITION  DOP 


MAIN 

2 WILL  BE  GREATER  THAN  NUMBER  LISTED', //»  5 

855  FORMATIlHl ,///, IX, 'DILUTION  OF  PRECISION  - ACCUMULATIVE  LATITUDE  0 I 

IISTRIBUTIQN' ,/, IX, 'PROBABILITY  THAT  FOUR  DIMENSIONAL  POSITION  OOP 
2WILL  BE  GREATER  THAN  NUMBER  LISTED',//! 

860  FORMAT ( IHl ,///, IX, ' DILUT ION  OF  PRECISION  PARAMETERS  - ACCUMULATIVE  j 

1 GLOBAL  01 STR IBJTION'  ,// ,8X, ' NUMBER* ,7X,  'VOOP' ,6X, 'HOOP' ,6X,  ! 

2 'MDOP' ,6X, 'TOOP'  ,6X, 'POOP* ,6X, 'GDOP' ,/ ) 

865  FORMATdH  , V X , F3  . I , 3 X , 6F  10. 4 ) 

870  F0RMAT(1H1,1X,'  MAXIMUM  AND  MINIMUM  NUMBERS  SEEN  AT  EACH  LATITUDE 
l£  LONGITUDE ' ) 

875  FORMATIIHO,/ ,38X,' MINIMUM' ,//,37X, 'LONGITUDE' ,//, 

1 32X,2  I • I' ,9X)  ,2C  2' , VX)  ,2(  ' 3'  ,VX)  ,/,22X,3(  • 5*  ,9X,'0' ,<JX)  , | 

2 '5' ,/ ,5X, 'LAT*,9X,8( ' O' ,9X! ,/ ) 

880  FORMATIIHO,/ ,38X,' MAXIMUM' ,// ,37X, 'LONGI TUOE  ' ,//, 

1 32X,2( ' 1*  ,9X) ,2 ('  2' ,9X1 ,2( *3'  ,9X) ,/,22X,3I ' 5* ,9X, *0* ,9X) , 

2 *5*,/,5X,  'LAT*,4X,8<'0'  ,9X1,/! 

885  FORMATdH  , 3X  , F9 .0 , 3X , 36  I 2 I 

890  FORMATdHl  ,90X,' NUMBER  OF  SATELL  I TE  S * , // , 5X  , *L  AT  ' ,3X  , 161  7 1 
895  FORMATdH  , 3X  , F9 .0 , 5X , 1 6F  7 . 2 ) 

900  FORMATdHO,/,lOX,' PROBABILITY  (IN  PERCENT)  OF  SEEING  EXACTLY  N SAT 
lELLITES'  ,/  I 

905  F0RMAT(1H0,/,10X, 'PROBABILITY  (IN  PERCENT)  OF  SEEING  N OR  MORE  SAT 
lELLITES* ,/ ) 

910  F0RMAT(1H0,/,1X, 'ON  A GLOBAL  BASIS  THE  PROBABILITY  (IN  PERCENT)  TH 
lAT  EXACTLY  N SATELLITES  WILL  BE  SEEN') 

915  FORMAT! IHO,/, IX, 'ON  A GLOBAL  BASIS  THE  PROBABILITY  (IN  PERCENT)  TH 
lAT  N OR  MORE  SATELLITES  WILL  BE  SEEN') 

920  CONTINUE 
END 


I 
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ORBINI 

SUBROUTINE  QRBINim 
C0MM0N/0RBIS/P(37,25I /CON/C( 18  J 
PU  ,2)  =»P(  1 .2  »*C(  1 ) 

P(I,3)=P(I  ,3)*C(l) 

P(I,<*)*P(I  1) 

P( 1 ,5»  = P( I t5  )*C( U 
P(  1 ,9»=S  INI  PU  ,<►)) 

PI  I ,10)=C0S(  PI  !,<»)) 

PI  I,m  = SINIPII,3»  I 
PI  I ,12 l=COSI PI  1,3) ) 

PIl  ,21)  = PI  I,22)*CI6) 

SAA  = 22808.*CI«>) 

PI  1 ,23)=IS AA/PI 1 ,21 ) )♦♦!  .5 

PI  I ,6I=PI I ,21)*I l.-PI I ,1)*PI  I,  1)  I 

PII,7)»S0RTI  CI12)/PII,6I) 

PII,8)>Cll2)/IPII,b)'*PII,6)) 

RETURN 

END 


T 


ORBIT 

SUBROUTINE  ORBI T ( I , T , PER , R , VE L, AC ) 

CGMMON/URBI S/OP(3/,2S)/CON/Cl  IB) 

DIMENSION  ROI.VELO),  AC  (3), TRAN  SO. 31,0(31 ,0EL(3»,0C(  3)  ,U(36I  , 

1 E(  36)  ,B  IGT(  36),  LILT(  36)  ,V(  36)  ,CV(  3o),  SV061  ,PER(361 
REAL  LILT 

IF(T.GT.O.  ) GO  TO  10 
Vm=0P(  1,5) 

SINE=(S0RT  (1  .-0P( I , I )**2)*S IN(0P( I ,5) ) )/ ( l.+0P( I , 1 ) *COS ( OP ( I , 5 ) ) ) 
RAO=(OP( I, 22)*(1 .-0P( I , 1 )**2) ) /( l.*0P( I , 1 ) *C OS ( OP ( I , 5 ) ) ) 

E(  I)=ARSIN(SINE) 

IF(0P( I, 1) .NE.O. ) GC  TO  15 
IFI V( I ) .LT . (C( 3)/2. ) ) E(I)»E(1) 

IF((V<I).GE.(C(3)/2.)).AND. 

1 (V(  I ) .LE.  1 3.*C<  3) /2.  ) ) ) E(  I )i=C(3)-E(I  ) 

lF(V(I).GT.(3.*C(3)/2.n  t(I)*C(^)tE(I) 

GO  TO  20 
15  CONTINUE 

IF(RAD.LT. OP( 1,22) .AND. V( I).L£.C(3))  E(I)=E(I) 

IF(RA0.GT.0P(I,22))  Em*CO)-E(I) 

IFIRAD.LT.UPI  I ,22  ) .AND.Vl  I I.GT  .C  (3)  ) E (I  ) =C  (<))♦£(  I ) 

20  CONTINUE 

BIGT(  n=  (PERI  I )*(E(  I )-0P(  I , 1 ) ♦SINE  ) )/(  2'».*C(  <i)  ) 

LILT(I)=BIGT(I) 

GO  TO  25 
10  DEL=.l 

LILT(1)  = BIGT(I)«-T 
E( I)=0. 

35  E( I )=E( I )FDEL 

Y=LILT  ( I l-(PER(I  )♦(  E(  I )-OP(  I.  1 )+SIN(  E(  I)  ))  ) / (2^.*C(^n 

IF(ABS(Y).LE.. 00001)  GO  TO  30 

IF(Y.GT.O.  ) GO  TO  35 

E( I)=E( I )-0EL 

0EL=DEL/10. 

GO  TO  35 
30  CONTINUE 

IF(E(I).GE.C(A))  E(  I)=E( I )-C(4l 

SINV=(  SORT (I .-OP( I , 1 )^^2 )*SIN(E( I ) ) )/( l.-OP(  I,  l)*COS(E ( I) ) ) 

V( I )=ARSIN(SINV) 

RP=OP( I ,22 )♦( 1 .-OP( I , 1 )^COS( E(  I ) ) ) 

P=OP( I ,22) ♦( l.-OP( I , I (♦♦2) 

IF(OP( I, 1) .NE.O. ) GO  TO  40 
IF(E(I).LT.(C(3)/2.)  ) V(I)=V(I) 

IF((E(1).GE.(C(3)/2.)>.AND. 

1 (E( I ) .L E. ( 3.^C(  3)/2. ) ) ) V(I )*C(3)-V(I  I 

IF(E( I ) .GT . ( 3.^C(3 )22. ) ) V { I)  = C { 4 ) ♦ V ( I ) 

GO  TO  25 
40  CONTINUE 

IF(RP.LT.P.AND.E(I).LE.C(4))  vm  = V(I) 

IF(RP.GT.P)  V( I )»C( 3)-V(  I ) 

IF(RP.LT.P  .AND.E(  I)  .GT.C(4)  ) V ( I ) =>C  ( 4 ) ♦V  ( I ) 

25  CONTINUE 

CV( I )=COS( V(  I ) ) 

SV( I )-SI N( V( I ) ) 

U(I)«=V(I)*OP(I,2) 

0P(I.13)=cC0S(U(I)) 

0P(I,14)==SIN(U(I  )) 

SPH»OP( I ,14) ♦0P( 1,9) 

CPH-S0RT(1 .tSPH^SPH) 


I 
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ORBIT 

OP(I ,15) =ARS  IN(SPH)*C (2) 

CPU  ,la)=U(  1 1 

Fl  = l.>OP(I  ,n*cvci) 

Om=OP(  I,6)/F1 

0(2)=0. 

Q(3)=0. 

OEL( l)=UP( 1 , 1)*QP( I , /)*SV(  I ) 
gEL(2)=0P< I ,7)*F1 
QEL(3)*0. 

QCm=-OPI  I,8)<'F1*FI 
QC(2) »0. 

0C(3)=0. 

CALL  TRMATX  (TRANS, I) 

CALL  MATMUL  (TRANS, Q,R) 

CALL  MATMUL  ( T RANS , 0 EL , V EL ) 

CALL  MATMUL  ( TRANS , QC , AC ) 

0P(1,2A)  = (T-  AINT(T«^.5)  )*C(^I  > 

0P( 1 ,16) =ATAN2(R (2 ) ,R( 1 ) )-0P( I ,2A) 
OP(I ,16)=0P(  I,16)*C(2) 

RETURN 

END 


! ^ 


i 

t 
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TRMATX 

SUBROUTINE  TRMATX(TR,N 
DIMENSION  TR(3,3) 

COMMON/ORBIS/OP( 37,25 )/CON/C( 18) 

TRU,l)*OP(I  ,12)*OP(  I,13)-0P(  I,m*OP{  1,  10I*OP(I,  U) 
TR II ,2)=OP( I ,12)*0P< I , 14)*0P( I , li)*OP( I, 10)*0P< I , 13) 
TR(l,3)*OPn,ll)*OP(  1.9) 

TR(2,1)=  OP(  I,ll)*OPn,13)+OP(  I,12)*0P(I,10)*0P(  I,U) 
TR(2,2)=-OP(  I,  11  )*OP(  I,  i^)«-aP(  I,  12  )*OP(  I , iO)*OP(  I ,13  ) 
TR(2,3)=-OP(  I, 12)*OP(I ,9) 

TR(  3,  l)*OP(  I ,9)*0P(  1 , 1<») 

TR(3,2)=  OP(  I,9)*0P(  1,13) 

TR(3,3)=OP( I ,10) 

TR(1,2)=-1.*TR( 1,2) 

RETURN 

END 


►1ATMUL 


SUBROUTINE  MAT MUL ( T , V , 0) 
DIMENSION  T(3,3> .VI3I ,0(3) 
DO  10  1*1,3 
0(1 )=0. 

DO  10  J=l,3 

10  0(I)*0(I)^T( I,J)*V(J) 
RETURN 
END 
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PQINT 

SUBROUTINE  POI NT ( AL D , AL A , T 1 M , VEC  ) 
CQMMUN/CON/C (181 
DIMENSION  VECI3) 

EW=ALO*C(l )*C(4)*TIM 
SN=ALA*C(1  ) 

VEC(  1 1 =C ( I OJ  *COS ( SN»  *COS ( E W» 
VEC(2)=C(lO)*COS(SN)*SIN(EW) 
VEC(3»=C(lO)*SIN(SN) 

RETURN 

END 


I 


i 

f 
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TAAT 


SUBROUTINE  TAAT  I MAX , MXX , MATR I X » 
DIMENSION  MATRIX(58905,4I 
DO  10  1=1,  58905 
00  10  II =l ,4 
10  MATRIX! I ,I I )=0 

IF(MAX.LT.3»  GO  TO  30 

MXX=(MAX-2»4(MAX-1 J*MAX/6 

MMMxMAX- I 

MMN=MAX 

MMO=MAX+l 

NA  = 0 

DO  25  K=2,MMM 
K0=K+1 

DO  20  L=KO,MMN 
KT=  L+1 

DO  15  M=KT,MMO 
NA=NAfl 

MATRIX(NA,2»  =K 
MATRIX(NA,3)=L 
MATRIX(NA,4>=M 
15  CONTINUE 
20  CONTINUE 
25  CONTINUE 
30  RETURN 
END 


1 

I 
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VOLUME 

SUBROUTINE  VOLUME  ( J V EC , I OSAT , VOL ) 

DIMENSION  UVECI37,3) , IDSAT  (4)  ,ONE(3)  ,TW0(3)  , THREEO)  .FOUR!  3) 
DIMENSION  TMF(3» ,TMT(3) ,0MT(3> .CROSS (3) 

KA=IDSAT(l  ) 

KB=IOSAT(2) 

KC=IOSAT (3  > 

KD  = IDSAT(<»  ) 

DO  10  N=l,  3 
ONE(N»=UVEC( KA.NI 
TWO(N)=UVEC(KB,N) 

THREE(N1=UVEC(KC.N> 

10  FOURIN)=UVEC(KO.N) 

CALL  VECTOR  ( TWO , 2 . F OUR , TMF I 
CALL  VECTOR  ( THR EE , 2 , TWO, TMT ) 

CALL  VECTOR  (ONE ,2 .TWO.OMT  ) 

CALL  VECTOR  ( TMT , 3 , TMF , C ROSS  I 
V0L=A8S ( DOT ( OMT, CROSS)) 

RETURN 

END 


J 

I 

i 


( 

) 

I 


I 
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COVNAV 

SUBROUTINE  COVNAV  ( G , I D , NA T , S I G I 

DIMENSION  I0(3/).S(^),SIGI6),G(37,4), IP  I VOT ( 4 ) , 1 NDE S ( 4, 2 ) 
REAL^o  TRA('»,^»tBBB(l,n  , DETER  M,  FPTMAX 
DATA  FPTMA X/Z7FFFFFFEFFFFFFFF/ 

DO  20  1=1, A 
DO  15  J=l, I 
TRAI I , J)=0. 

00  10  K=1,NAT 
L=ID( K) 

TRA(I,J)=TRAri,J)+G(L,I)*G(L,J) 

10  CONTINUE 

TRA(J,n=TRA(  I,J» 

15  CONTINUE 

TRAI  I , I ) =TRA  ( I ,I  )H  .E-12 
20  CONTINUE 

CALL  DMATRI  I TR A , A , A , BBB , 1 , 0 , I P I VOT , I NOE S , I S I NG , DETE RM ) 

IF! ISING.NE.OI  STOP 
IFIDETERM.EO. FPTMAX)  GO  TO  25 
SIGl I )=DSORT (TRA(3,3  ) ) 

S IGI2 )=DSORT (TRAI I , 1) fTRA(2,2  ) ) 

S IGI 3)=0MAXl (OSyRTITRA I 1 ,1 ) ) , DSOR T I TRA 1 2 , 2 ) ) ) 

SIGIA)=D50RT ITRAIA.A) ) 

S IGI  5 )=D  SORT  (TRA  II  ,n+TRAI2,2)»TRAI3,3n 
SIGIb)=0SQRTITRAIl,n+TRAI2,2»fTRAI3,3)frTRAIA,A)) 

RETURN 

25  WRITEI6.30) 

30  FORMATIIH  , / , 3X,  ' D ET E RM I N ATE  REACHED  MAXIMUM  VALUE*) 

STOP 

END 


! 


I 


BLK  DATA 

BLOCK  DATA 
CCMMON/CON/CCIB) 

DATA  C /. 01 745329252, 57. 295 7795 13, 3. U1 59 26536, 6. 2 831 85 307 18, 
11 .57  0 79630, 6 0 7b.  11  6,  0. ,0. , 0. , 2 . 09 26 1 43 504E ♦ 0 7 , 

27.2921 1585E- 05, 1.40 76380E ♦16,365.2563835  ,92.91 
3E ♦06,0. 016  72  72,23.44436,-  77.7303,5280. / 

END 
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MICROCOPY  RESOLUTION  TEST  CHART 
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^ 1 


I VECTOR 

SUBROUTINE  VECTOR! VI , 1 ,V2*V3) 
DIMENSION  V1(3).V2(3),V3(3) 

GO  TO  (10, 15,20) . I 
10  v3m<=vi(i)«^v2m 
V3(2)»V1(2)«-V2(2) 

, V3(3)*Vi(3)«-V2(3) 

I RETURN 

15  V3(1)>V1(1)-V2(1) 
V3(2)>>V1(2)-V2(2) 
V3(3)*V1(3)-V2(3) 

RETURN 

'l  20  V3(1)*VI(2)*V2(3)-V1(3I*V2(2I 

V3(2)>V1(3)*V2(1)-V1( 1)*V2(3) 
V3(3)*V1( 1)*V2(2)-V1(2)*V2(1) 
RETURN 
END 


\ 


I 


i 


r 


FUNCTION  OOT(Vl.V2l 

; OIMENSIONVU3I  tV2(3l 

I DaT«Vl(l}*V2m«^Vl(2)*V2(2)«-Vl(3l*V2(3l 

I RETURN 

I END 

I 

I 

4 
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UNiVEC 

SUBROUTINE  UNIVEC  (VtUV) 

DIMENSION  V(3)tUV(3l 
OENOM  » SORT (DOT (Vt VI > 

UVUJ  » V(  II /OENOM 
UVI21  * VI2I/DENOM  1 i 

UV(3I  » V(3)/OENOM  i | 

RETURN  II 


^ I e iktr% 


SUBROUTINE  TALL  (NAX.NXX.MATRIX} 
DIMENSION  NATRIX(58905t4l 
DO  10  I«lt5»405 
DO  10  11*1.4 
10  MATRIxn*lI)-0 

IFIMAX.lt. 31  GO  TO  35 

MXX-(MAX-3)*(MAX-2l*fMAX-l)*MAX/24 

KK«MAX-3 

LL-MAX-2 

MM«MAX-1 

NN«MAX 

NA*0 

DO  30  K>1.KK 
K0*K^1 

DO  25  L»KO.LL 
KT*L+1 

DO  20  M«KT.MM 
KP»M*1 

DO  15  N*KP.NN 
NA«NA4-1 

MATRIX(NA,1)-K 
MATRIXCNA.21-L 
MATRIX(NA.3)-M 
MATR1X(NA.4)*N 
15  CONTINUE 
20  CONTINUE 
25  CONTINUE 
30  CONTINUE 
35  RETURN 
END 
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